{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating Confidence Intervals for Machine Learning Classifiers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Author: Sebastian Raschka\n",
      "\n",
      "Last updated: 2022-04-24\n",
      "\n",
      "Python implementation: CPython\n",
      "Python version       : 3.9.7\n",
      "IPython version      : 8.0.1\n",
      "\n",
      "numpy     : 1.22.1\n",
      "mlxtend   : 0.19.0\n",
      "matplotlib: 3.5.1\n",
      "sklearn   : 1.0.2\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -a 'Sebastian Raschka' -u -d -v -p numpy,mlxtend,matplotlib,sklearn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- This notebook provides an outline of different methods for creating confidence intervals for machine learning models. Note that these methods also apply to deep learning. \n",
    "- The text in this notebook is purposefully short so that we can focus on the technical execution without getting bogged down in details; there are many links to all the relevant conceptual explanations throughout this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence Intervals in a Nutshell"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- In a nutshell, what is a 95% confidence interval anyway? Imagine that you have a statistic (a measurement of some sort) that you calculated from a sample drawn from a population. Our goal is to estimate a population parameter with this statistic, for example, we could estimate the population mean with the sample mean. Most of the time, the estimated value and the true value are not the same, and we can use the confidence interval to quantify the uncertainty of that estimate.\n",
    "\n",
    "- But how do we interpret confidence interval? A 95% confidence interval is essentially a method that computes an upper and lower bound around that estimate. The true parameter is either insider or outside these bounds. \n",
    "\n",
    "\n",
    "\n",
    "- In practice, this is not feasible, but let's assume for a moment that we had access to the population. If we drew a very large number of samples from the distribution and applied our confidence interval method to these samples, 95% of the confidence intervals would contain the true value.\n",
    "\n",
    "- In a machine learning context, what we are usually interested in is the performance of our model. So, here the population parameter we want to estimate could be the generalization accuracy of our model. Then, the test set accuracy represents our estimated generalization accuracy. The 95% confidence interval gives us an uncertainty measure of how accurate this estimate is.\n",
    "\n",
    "- As a side-note, we can say that the difference of two measurements is statistically significant if confidence intervals *do not* overlap. However we *cannot* say that results are *not* statistically significant if confidence intervals overlap. (The [Error bars](https://www.nature.com/articles/nmeth.2659) article in the Points of Significance series illustrates this nicely.) If we want to check whether the difference is not statistical significant, we would have to take a look at the distribution of the differences and check whether its confidence interval contains 0 or not. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining a Dataset and Model for Hands-on Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- We just use the [Iris dataset](https://archive.ics.uci.edu/ml/datasets/iris) and a [Decision tree classifier](https://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeClassifier.html) for the sake of simplicity.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mlxtend.data import iris_data\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "\n",
    "X, y = iris_data()\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X, y, test_size=0.15, random_state=123, stratify=y\n",
    ")\n",
    "\n",
    "clf = DecisionTreeClassifier(random_state=123)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Normal Approximation Interval Based on the Test Set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- This is maybe the easiest and classic way for creating confidence intervals.\n",
    "\n",
    "- Using this method, we compute the confidence interval from a single training-test split. This is particularly attractive in deep learning where model training is expensive.\n",
    "\n",
    "- It's also attractive (usually in a deep learning context) when we are interested in a very particular model (vs. models fit on different training folds like in k-fold cross-validation).\n",
    "\n",
    "- You can find a description of these method i section *1.7 Confidence Intervals via Normal Approximation* of my \"[Model Evaluation, Model Selection, and Algorithm Selection in Machine Learning](https://arxiv.org/abs/1811.12808)\" article."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- In a nutshell, the equation for computing the confidence interval for an estimated parameter (let's say the the sample mean $\\bar{x}$) assuming a normal distribution is computed as follows:\n",
    "\n",
    "$$\n",
    "\\bar{x} \\pm z \\times \\text{SE}\n",
    "$$\n",
    "\n",
    "- where\n",
    "    - $z$ is the $z$ value (the number of standard deviations that a value lies from the mean of a standard normal distribution)\n",
    "    - $SE$ is the standard error of the estimated parameter (here: sample mean\n",
    "    \n",
    "    \n",
    "- In our case, the sample mean $\\bar{x}$ is test set accuracy $\\text{ACC}_{\\text{test}}$, a proportion of success (in the context of a [Binomial proportion confidence interval](https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval)). \n",
    "- The standard error, under a normal approximation can be computed as \n",
    "\n",
    "\n",
    "$$\n",
    "\\text{SE} = \\sqrt{ \\frac{1}{n} \\text{ACC}_{\\text{test}}\\left(1- \\text{ACC}_{\\text{test}}\\right)},\n",
    "$$\n",
    "\n",
    "where $n$ is the test set size.\n",
    "\n",
    "- So, plugging the SE back into the formula above, we get\n",
    "\n",
    "$$\n",
    "\\text{ACC}_{\\text{test}} \\pm z \\sqrt{\\frac{1}{n} \\text{ACC}_{\\text{test}}\\left(1- \\text{ACC}_{\\text{test}}\\right)}.\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Let's start with computing the z-value, which can obtain from `scipy.stats.norm.ppf` rather than looking it up from a $z$ table in one of our old statistics textbook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.959963984540054\n"
     ]
    }
   ],
   "source": [
    "import scipy.stats\n",
    "\n",
    "confidence = 0.95  # Change to your desired confidence level\n",
    "z_value = scipy.stats.norm.ppf((1 + confidence) / 2.0)\n",
    "print(z_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Next, let's compute the test accuracy of the classifier, and plug in the values into the formula above; the Python code for this is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.873179017733963 1.0398644605269067\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "clf.fit(X_train, y_train)\n",
    "\n",
    "acc_test = clf.score(X_test, y_test)\n",
    "ci_length = z_value * np.sqrt((acc_test * (1 - acc_test)) / y_test.shape[0])\n",
    "\n",
    "ci_lower = acc_test - ci_length\n",
    "ci_upper = acc_test + ci_length\n",
    "\n",
    "print(ci_lower, ci_upper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Let's visualize the confidence interval using the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAADQCAYAAAD4dzNkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUC0lEQVR4nO3dfbRkVX3m8e8jYEBAGpT0QjRgJLZxwGkIYUAMNurEyDICxiTOMCMQUcjEqBmRRYKjLCejRKNOFIWIQcVlDEpAAVcmJLwmKCCv3byqEdSAA0FsQAYB8Td/nH2xuNyX7svtW7Xx+1mr1q06dV5+tet0PbX3OV0nVYUkSerLk8ZdgCRJWn8GuCRJHTLAJUnqkAEuSVKHDHBJkjq08bgLkGazbNmy2mmnncZdxmPcd999bL755uMuY0aTWtuk1gXWthCTWhdMbm1XXHHFnVW17WKu0wDXxFq+fDmXX375uMt4jAsuuIBVq1aNu4wZTWptk1oXWNtCTGpdMLm1Jfn2Yq/TIXRJkjpkgEuS1CEDXJKkDhngkiR1yACXJKlDBrgkSR0ywCVJ6pABLklShwxwSZI6ZIBLktQhA1ySpA4Z4JIkdcgAlySpQwa4JEkdMsAlSeqQAS5JUocMcEmSOmSAS5LUIQNckqQOGeCSJHXIAJckqUMGuCRJHTLAJUnqkAEuSVKHDHBJkjpkgEuS1CEDXJKkDhngkiR1yACXJKlDBrgkSR0ywCVJ6pABLklShwxwSZI6ZIBLktQhA1ySpA4Z4JIkdcgAlySpQwa4JEkdMsAlSeqQAS5JUocMcEmSOmSAS5LUIQNckqQOGeCSJHXIAJckqUMGuCRJHTLAJUnqkAEuSVKHDHBJkjpkgEuS1CEDXJKkDhngkiR1yACXJKlDBrgkSR0ywCVJ6pABLklShwxwSZI6ZIBLktQhA1ySpA4Z4JIkdcgAlySpQwa4JEkdMsAlSeqQAS5JUocMcEmSOmSAS5LUIQNckqQOGeCSJHXIAJckqUMGuCRJHTLAJUnqkAEuSVKHDHBJkjpkgEuS1CEDXJKkDhngkiR1yACXJKlD8wZ4kkrygZHHRyY5doNW9dgaLkiy+1Juc30keXeSly3CepYl+W8jj5+R5LTHu962rk8kef488xww3zyLVMshSY7f0NuRNB5fvOpW9j7uPJ599JfZ+7jz+OJVt467pCekdemBPwC8OsnTF7KBJBsvZLlxSbLR+i5TVe+sqn9chM0vAx4J8Kq6rapeswjrpaoOq6rr55ntAGC9Ary391fShvXFq27lj09fw61r76eAW9fezx+fvsYQ3wDW5cP3x8DHgT8Cjhl9IskOwMnAtsC/AYdW1XeSfAq4C9gVuDLJ04D7gecBOwCHAgcDewGXVtUhbX0nAL8KbAacVlXvmquwJO8EfrPN/xXg8KqqJBcAVwN7AE8Ffq+qLmsjB88BtgeeBbyvqk5Ksgp4F/A9YGWS3YATgN3b6//vVXV+ki8Bf1tVpyQ5HNinqg5qr/fsqjotyS3AXwP7ApsAbwTeC+wEvL+qTkyyBfAlYOs2zzuq6kvAccBzklwN/APw0bbenZNsOktNhwCvAp7SXtsZVXXUDG11AXBkVV2e5IfAXwCvbO/L/m3ZVwEvTvIO4Lfaoh9leH//H/CGqrpx2vt7dZIDgZVVtbZt65vA3q393wE8Gfg+cFBV3T7Xe6onlt/9y68CsHbt/Zxw01fHXM3MrG39zVXXVd9Zy4MP/+RR0+5/6GGOOm01n7vsO2Ot7YlmXXtPHwVWJ3nftOnHA6dU1aeT/B7wYYZeHMBzgZdV1cPtA39r4CUMIXEWwwf8YcDXkqysqquBY6rqrtYLPjfJC6pq9Rx1HV9V7wZI8hmGQDqrPbd5Vb0wyT4MXzJ2btNfAOwJbA5cleTLbfoewM5VdXOStwFU1S5Jngeck+S5DGF8cZKbgbe19czku1W1V5IPAZ9qr3VT4DrgROBHwIFVdU8b2bgkyZnA0a2Gle017Tiyzj+YpSaAlQxh+gBwU5KPVNV352i3zYFLquqY9p6+oar+tNVwdlWd1rZ/LnBEVX0jyX8APsbwHsKj398nAQcCn2zz3VJVtyf5Z2DP9qXqMOCo1m6zSvJGhnZm+fLlc80qaQJND+/5pmvh1inAW9CcAryZocc2ZS/g1e3+Z4DRgP9CVT088vis9kG+Bri9qtYAJLkO2JGhx/w77QN8Y2A7huHcuQJ83yRHMfQ+t2EIyKkA/1yr/aIkT02yrE3/UlXdD9yf5HyG4F4LXFZVN7d5XgR8pC1/Y5JvA8+tqtWt138+QwDfNUtdZ7a/a4Atqupe4N4kP2p13Ae8p325+AnDiMB8aTVjTe25c6vqboAk1zOMcswV4A8CZ7f7VwD/cfoMbZTghcAXkkxN/rmRWUbf31OBdwKfBF7bHgM8Ezg1yXYMvfCbmUdVfZxhxIcVK1bUfPNrsp16+F4AXHDBBaxatdeYq5mZta2/uera+7jzuHXt/Y+Zvv2yzR7ZHzakSW2zzx+x+Otcn7PQ/zfweobe22xGP3Dvm/bcA+3vT0buTz3eOMmzgSOBl1bVC4AvM/RaZ9SGlD8GvKaqdgFOmjb/9A//mmf6aL1hdrswDAc/Y4555nytwEEMw9K/0nrbtzPHa12Hmka38TDzfzF7qKqmXvds8z8JWFtVK0duvzzy/Gh7fRXYKcm2DCMwp7fpH2EYJdkFOJz5X6Okzr395SvYbJNHn0q02SYb8faXrxhTRU9c6xzgrbf5eYYQn/IVhh4XDKH0z4+jlqcyhMLdSZYDr5hn/qkwuLP1Fqef7PW7AEleBNw91UMF9k+yaTsuvwr42gzrvojh9dCGqX+BYWh6j1bXrsCR7UvHQmwF3FFVDyXZl6HHDHAvsOUsy8xY0wK3P5tHtl9V9wA3J/ntts0k+fczLdS+DJwBfBC4oaq+357aCpg6c+XgRa5V0gQ6YNftee+rd2H7ZZsRhp73e1+9Cwfsuv24S3vCWd8ziD8AvGnk8ZuBk5O8nXYS20ILqaprklzFMAz+LeDieeZfm+QkhmHqW3hsEP8gyVdoJ7GNTL+MoXf/C8D/rKrbRo4lT/kYcGIb7v8xcEibfhLDiXq3tePkJyd5Cevvs8BZSS5nOHRwY3tN309ycZJrgb9jOPdg1pqq6oGR4e3F8DfASUnezPCF6CDghHZS2ybt+WtmWfZUhvfgkJFpxzIMwd8KXAIs9AuPpI4csOv2BvYSyE9HUp84Rs+4njb9WOCHVfXn46hL62fFihV1002LPcjw+A3H2FaNu4wZTWptk1oXWNtCTGpdMLm1Jbmiqhb190z8JTZJkjr0hPwRjqpaNcv0Y5e2EkmSNgx74JIkdcgAlySpQwa4JEkdMsAlSeqQAS5JUocMcEmSOmSAS5LUIQNckqQOGeCSJHXIAJckqUMGuCRJHTLAJUnqkAEuSVKHDHBJkjpkgEuS1CEDXJKkDhngkiR1yACXJKlDBrgkSR0ywCVJ6pABLklShwxwSZI6ZIBLktQhA1ySpA4Z4JIkdcgAlySpQwa4JEkdMsAlSeqQAS5JUocMcEmSOmSAS5LUIQNckqQOGeCSJHXIAJckqUMGuCRJHTLAJUnqkAEuSVKHDHBJkjpkgEuS1CEDXJKkDhngkiR1yACXJKlDBrgkSR0ywCVJ6pABLklShwxwSZI6ZIBLktQhA1ySpA4Z4JIkdcgAlySpQwa4JEkdMsAlSeqQAS5JUocMcEmSOmSAS5LUIQNckqQOGeCSJHXIAJckqUMGuCRJHTLAJUnqkAEuSVKHDHBJkjpkgEuS1CEDXJKkDhngkiR1yACXJKlDBrgkSR0ywCVJ6pABLklShwxwSZI6ZIBLktQhA1ySpA4Z4JIkdcgAlySpQ6mqcdcgzSjJvcBN465jBk8H7hx3EbOY1NomtS6wtoWY1LpgcmtbUVVbLuYKN17MlUmL7Kaq2n3cRUyX5PJJrAsmt7ZJrQusbSEmtS6Y3NqSXL7Y63QIXZKkDhngkiR1yADXJPv4uAuYxaTWBZNb26TWBda2EJNaF0xubYtelyexSZLUIXvgkiR1yACXJKlDBriWRJLfSHJTkm8mOXqG57dKclaSa5Jcl+TQ+ZZNsk2Sf0jyjfZ366WsLcmzkpyf5IY2/S0jyxyb5NYkV7fbfktVV3vuliRr2rYvH5k+7jZbMdImVye5J8lb23NL0WZbJzkjyeoklyXZeb5ll7DNZqxtAvazudps3PvZbG22ofezk5PckeTaWZ5Pkg+3ulcn2W2+17SgNqsqb9426A3YCPgX4BeBJwPXAM+fNs+fAH/W7m8L3NXmnXVZ4H3A0e3+0VPLL2Ft2wG7telbAl8fqe1Y4MhxtFl7fAvw9BnWO9Y2m2E9/xfYYQnb7P3Au9r95wHnzrfsErbZbLWNez+bsa4J2c9mrW1D7WdtHfsAuwHXzvL8fsDfAQH2BC7dEPuZPXAthT2Ab1bVt6rqQeBvgP2nzVPAlkkCbMHwgf/jeZbdH/h0u/9p4IClrK2qvldVVwJU1b3ADcD2C6hhUeuaZ71jbbNp87wU+Jeq+vYCalhoXc8HzgWoqhuBHZMsn2fZpWqzGWubgP1stjaby1jbbNo8i72fUVUXMezTs9kfOKUGlwDLkmzHIu9nBriWwvbAd0ce/yuP/QA6Hvhl4DZgDfCWqvrJPMsur6rvAbS/P7/EtT0iyY7ArsClI5Pf1IbPTl7AEOLjrauAc5JckeSNI8tMTJsBrwU+N23ahm6za4BXAyTZA9gBeOY8yy5Vm81W2yPGtJ/NVde497N524zF38/WxWy1L+p+ZoBrKWSGadP//+LLgauBZwArgeOTPHUdlx1XbcMKki2AvwXeWlX3tMknAM9p838P+MAS17V3Ve0GvAL4gyT7rOf2N2RtJHky8CrgCyPLLEWbHQdsneRq4A+BqxhGBiZhP5uttmEF49vP5qpr3PvZfG22IfazdTFb7Yu6nxngWgr/Cjxr5PEzGXpmow4FTm9DTt8EbmY4pjXXsre3YSna3zuWuDaSbMLwofrZqjp9aoGqur2qHm69zpMYhs6WrK6quq39vQM4Y2T7Y2+z5hXAlVV1+9SEpWizqrqnqg6tqpXA6xiOz988z7JL0mZz1DbW/Wyuusa9n81VW7Mh9rPHU/ui7mcGuJbC14BfSvLs9o34tcCZ0+b5DsOxKtoxrBXAt+ZZ9kzg4Hb/YOBLS1lbO777V8ANVfXB0QWm/iE2BwIznq26geraPMmWbfrmwK+PbH+sbTby/H9i2rDmUrRZkmXtOYDDgItab3bs+9lstY17P5ujrrHvZ3O8n1M2xH62Ls4EXpfBnsDdbVh8cfezuc5w8+ZtsW4MZ2V+neEMzGPatCOAI9r9ZwDnMBwvvRb4L3Mt26Y/jeEElm+0v9ssZW3AixiGv1YzDBdfDezXnvtMm391+4e53RLW9YsMxwavAa6bpDZrzz0F+D6w1bR1LkWb7dVe+43A6cDWE7SfzVjbBOxns9U1CfvZXO/nhtzPPscw/P4QQ6/69dPqCvDRVvcaYPcNsZ/5U6qSJHXIIXRJkjpkgEuS1CEDXJKkDhngkiR1yACXJKlDBrikxyXJwxmu6nRtki8kecrjWNenkrym3f9EkufPMe+qJC8ceXxEktctdNtSbwxwSY/X/VW1sqp2Bh5k+P+wj0iy0UJWWlWHVdX1c8yyCngkwKvqxKo6ZSHbGpckG4+7BvXLAJe0mP4J2Kn1js9P8tfAmiQbJXl/kq+1i0gcDo9cN/n4JNcn+TIjF3BIckGS3dv930hyZYbri5+b4aIeRwB/1Hr/v5bhOs9HtvlXJrmkbeuMtAtWtHX+WYZrR389ya9NfwFJtmjbuDLDta73H3nudW2d1yT5TJu2vG3jmnZ7YZIdM3Kt6CRHJjl2pIb3JLkQeEuS30xyaZKrkvxj++W6qTo+2WpYneS3krw+yYdG1vuGJI/6dTb97PDbn6RF0XqTrwD+T5u0B7BzVd2c4UpVd1fVryb5OeDiJOcwXFlrBbALsBy4Hjh52nq3ZfjN6n3aurapqruSnAj8sKr+vM330pHFTgH+sKouTPJu4F3AW9tzG1fVHkn2a9NfNu2l/Ag4sIafC306cEmSMxkuXXkMwwU87kyyTZv/w8CFVXVgG23YApjvClfLqurFre6tgT2rqpIcBhwFvA34H63NdhmZ70FgdZKjquohht+cP3yebekJygCX9HhtluFqUDD0wP+KYWj7sqqaurDErwMvmDq+DWwF/BKwD/C5qnoYuC3JeTOsf0+G37ieuoDGXNdhJslWDAF5YZv0aR59Naqpi4FcAew40yqA92S4stbUJW2XAy8BTquqO6fV8RKGC2nQXsfdmf8SlaeO3H8mcGr7je4n89OLcbyM4beyaev+QXt95wGvTHIDsElVrZlnW3qCMsAlPV7313A1qEckAbhvdBJDj/jvp823H/NfTjHrMM/6eKD9fZiZPwMPYriq1a9U1UNJbgE2Xc86fsyjD1FuOu350bb5CPDBqjozySrg2DZ9tu19AvgTht///uQ61qMnII+BS1oKfw/8fobLYpLkuRmuYHUR8Np2jHw7YN8Zlv0q8OIkz27LTg1d3wtsOX3mqrob+MHI8e3/Clw4fb45bAXc0cJ7X2CHNv1c4HeSPG1aHecCv9+mbZThuue3Az+f5GntkMEr59nere3+wSPTzwHeNPVgqldfVZcyXJLyPzPtSlv62WKAS1oKn2A4vn1lO7nrLxl6v2cwXH1pDXACMwRtVf0b8Ebg9CTX8NPh57OAA6dOYpu22MHA+5OsBlYC716PWj8L7J7kcobe+I2tjuuA/wVc2OqYOnnsLcC+SdYwDMv/u3Z8+t3ApcDZU+uYxbHAF5L8E3DnyPQ/BbbO8N/zruHRX24+D1w8Nayun01ejUySOpPkbOBDVXXuuGvR+NgDl6ROJFmW5OsM5x0Y3j/j7IFLktQhe+CSJHXIAJckqUMGuCRJHTLAJUnqkAEuSVKH/j8r7bkDTEHpzAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(7, 3))\n",
    "\n",
    "ax.errorbar(acc_test, 0, xerr=ci_length, fmt=\"o\")\n",
    "\n",
    "ax.set_xlim([0.8, 1.0])\n",
    "\n",
    "ax.set_yticks(np.arange(1))\n",
    "ax.set_yticklabels([\"Normal approximation interval\"])\n",
    "ax.set_xlabel(\"Prediction accuracy\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.grid(axis=\"x\")\n",
    "plt.savefig(\"matplotlib-figures/normal-approx.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Lastly, let's store our confidence interval in a Python dictionary so that we can retrieve it later when we compare it to other confidence intervals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = {\n",
    "    \"Method 1: Normal approximation\": {\n",
    "        \"Test accuracy\": acc_test,\n",
    "        \"Lower 95% CI\": ci_lower,\n",
    "        \"Upper 95% CI\": ci_upper,\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2) Bootstrapping Training Sets -- Setup Step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Confidence intervals are used to estimate unknown parameters. If we only have one measurement, like the accuracy from a single test set, we need to make assumptions about the distribution that this accuracy values comes from. E.g., we can assume it's from a normal distribution. \n",
    "\n",
    "- In an ideal world, we have access to the distribution that this sample comes from, so that we can e.g., look at the range of values that 95% of the samples accuracy values fall into. This is desirable but not practical since we don't have an infinite pool of test sets.   \n",
    "\n",
    "- Now, a workaround is bootstrapping, which estimates the sampling distribution. This is done by taking multiple samples with replacement from a single random sample.\n",
    "\n",
    "- The method we are using here is also often referred to as out-of-bag bootstrap where we are training the model on training folds and evaluate it on held-out data points from each round. For more detail, please see section 2, *Bootstrapping and Uncertainties* of \"[Model Evaluation, Model Selection, and Algorithm Selection in Machine Learning](https://arxiv.org/abs/1811.12808).\"\n",
    "\n",
    "- The equation is as follows:\n",
    "\n",
    "$$\n",
    "\\text{ACC}_{\\text{bootavg}}=\\frac{1}{b} \\sum_{j=1}^{b} \\text{ACC}_{\\text{boot}, j}\n",
    "$$\n",
    "\n",
    "- There are different flavors of bootstrapping. Below, we create multiple bootstrap samples for re-use in the upcoming sections.\n",
    "\n",
    "- Note that a minimum of bootstrap rounds of 200 is the typical recommendation as per the \"[Introduction to the Bootstrap](https://scholar.google.com/scholar_lookup?title=An%20introduction%20to%20the%20bootstrap&author=&publication_year=1993)\" book. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9463377985233019"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "rng = np.random.RandomState(seed=12345)\n",
    "idx = np.arange(y_train.shape[0])\n",
    "\n",
    "bootstrap_train_accuracies = []\n",
    "bootstrap_rounds = 200\n",
    "\n",
    "\n",
    "for i in range(bootstrap_rounds):\n",
    "\n",
    "    train_idx = rng.choice(idx, size=idx.shape[0], replace=True)\n",
    "    valid_idx = np.setdiff1d(idx, train_idx, assume_unique=False)\n",
    "\n",
    "    boot_train_X, boot_train_y = X_train[train_idx], y_train[train_idx]\n",
    "    boot_valid_X, boot_valid_y = X_train[valid_idx], y_train[valid_idx]\n",
    "\n",
    "    clf.fit(boot_train_X, boot_train_y)\n",
    "    acc = clf.score(boot_valid_X, boot_valid_y)\n",
    "    bootstrap_train_accuracies.append(acc)\n",
    "\n",
    "bootstrap_train_mean = np.mean(bootstrap_train_accuracies)\n",
    "bootstrap_train_mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- We can visualize the test accuracies from bootstrapping ($ \\text{ACC}_{\\text{boot}, j}$) along with their sample mean ($\\text{ACC}_{\\text{bootavg}}$) in a historgram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEYCAYAAABRMYxdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAb/klEQVR4nO3de5CdBZnn8e+TADYkkUBsMoFQkzBQrQ5i0HATd6oxODiKJmOJKzsjUdmNU64Ozq6s4Bar69YWTA21NS64apZxbFgGjVyWMOM6ZKLHy5oJJFwCGAKOhNhDb6IJARqSmMuzf5wXDOnu5KT7vOd0v+f7qeo6570//fh2+PleIzORJEmqkkntLkCSJKnZDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyjmh3AY2YPn16nnrqqe0uo1JefPFFpkyZ0u4yKsWeNp89bT572nz2tLnWrl37q8zsHut6JkTAmTlzJmvWrGl3GZVSq9Xo7e1tdxmVYk+bz542nz1tPnvaXBHxdDPW4ykqSZJUOQYcSZJUOQYcSZJUORPiGpzh7N69m/7+fnbu3NnuUlqqq6uL2bNnc+SRR7a7FEmSxq0JG3D6+/uZNm0ac+bMISLaXU5LZCZbt26lv7+fuXPntrscSZLGrQl7imrnzp3MmDGjY8INQEQwY8aMjjtqJUnS4So14ETEn0XEYxHxaETcFhFdEXF8RKyIiCeLz+PGsP5mljshdOLvLEnS4Sot4ETEScCfAvMz83RgMvAh4CpgZWaeBqwshiVJkpqm7FNURwBHR8QRwDHAM8BCoK+Y3gcsKrmG0kQEH/7wh18Z3rNnD93d3Vx88cVtrEqSJJUWcDLzn4HrgU3AAPBcZt4LzMzMgWKeAeCEsmoo25QpU3j00UfZsWMHACtWrOCkk05qc1WSJKm0u6iKa2sWAnOB7cC3I+KPD2P5JcASgO7ubmq12qumH3vssbzwwgvNKnfUFixYwO23386iRYu4+eabef/7389PfvITXnjhBV588UWuvPJKHnvsMfbu3cvVV1/Ne97zHp5++mmWLFnCSy+9BMD111/POeecw49+9COuvfZaZsyYwU9/+lPmzZvHTTfdNOS6m507dw7px+EaHBwc8zr0avb08P2ofze/2pG87ujgX8we+ugDe9p89rT57On4VOZt4hcCT2XmLwEi4k7gbcDmiJiVmQMRMQvYMtzCmbkUWArQ09OTB77nY/369UybNu2V4W+v+QW3r+1vWvEfeOtsLpl/8iHnu+yyy/jiF7/IJZdcwvr16/n4xz/Offfdx7Rp07j22mu56KKLuOWWW9i+fTtnn302733veznllFP43ve+R1dXF08++SSXXnopa9as4ZhjjmHdunU89thjnHjiiZx//vmsW7eOt7/97a/aZldXF2eeeeaYfj/fndJ89vTwfeVrq1j91DbOmXs81/SeN2S6PW0+e9p89nR8KjPgbALOjYhjgB3AAmAN8CKwGLiu+Ly7GRvrf3YHq5/a1oxVAXDuKTMamu+MM85g48aN3Hbbbbz73e9+1bR7772X5cuXc/311wP1Iy+bNm3ixBNP5JOf/CQPPfQQkydP5oknnnhlmbPPPpvZs2cDMG/ePDZu3Dgk4EiSpIMrLeBk5uqIuB14ANgDPEj9iMxUYFlEXE49BF3SjO3NPu5ozpl7fDNW9cr6GvW+972Pz3zmM9RqNbZu3frK+MzkjjvuoKen51Xzf+ELX2DmzJk8/PDD7Nu3j66urlemveY1r3nl++TJk9mzZ88YfgtJkjpTqU8yzszPA58/YPQu6kdzmuqS+Sc3dEqpDB/72Mc49thjedOb3vSq87AXXXQRN9xwAzfccAMRwYMPPsiZZ57Jc889x+zZs5k0aRJ9fX3s3bu3LXVLklRVE/ZJxuPJ7NmzueKKK4aMv+aaa9i9ezdnnHEGp59+Otdccw0An/jEJ+jr6+Pcc8/liSeeYMqUKa0uWZKkSpuw76IaDwYHB4eM6+3tfeVis6OPPpqvfe1rQ+Y57bTTWLdu3SvD11577ZBlAW688cbmFixJUofwCI4kSaocA44kSaqcCR1wMrPdJbRcJ/7OkiQdrgkbcLq6uti6dWtH/Qc/M9m6deurbiuXJElDTdiLjGfPnk1/fz+//OUv211KS3V1db3yIEBJkjS8CRtwjjzySObOndvuMiRJ0jg0YU9RSZIkjcSAI0mSKseAI0mSKseAI0mSKseAI0mSKseAI0mSKseAI0mSKseAI0mSKseAI0mSKseAI0mSKseAI0mSKqe0gBMRPRHx0H4/z0fEpyPi+IhYERFPFp/HlVWDJEnqTKUFnMzckJnzMnMe8FbgJeAu4CpgZWaeBqwshiVJkpqmVaeoFgD/lJlPAwuBvmJ8H7CoRTVIkqQO0aqA8yHgtuL7zMwcACg+T2hRDZIkqUNEZpa7gYijgGeA383MzRGxPTOn7zf92cwcch1ORCwBlgB0d3e/ddmyZaXW2WkGBweZOnVqu8uoFHt6+K5dvYMNz+6j57hJXH3O0UOm29Pms6fNZ0+b64ILLlibmfPHup4jmlHMIfwB8EBmbi6GN0fErMwciIhZwJbhFsrMpcBSgJ6enuzt7W1BqZ2jVqthT5vLnh6+r2xYBc9uY/r06fT2njdkuj1tPnvafPZ0fGrFKapL+c3pKYDlwOLi+2Lg7hbUIEmSOkipAScijgHeCdy53+jrgHdGxJPFtOvKrEGSJHWeUk9RZeZLwIwDxm2lfleVJElSKXySsSRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqpxSA05ETI+I2yPi8YhYHxHnRcTxEbEiIp4sPo8rswZJktR5jih5/V8CvpuZH4iIo4BjgM8BKzPzuoi4CrgK+GzJdUhSS9yzod0VHNy0dhcgtUhpR3Ai4rXA7wF/BZCZv87M7cBCoK+YrQ9YVFYNkiSpM5V5iuoU4JfAX0fEgxFxU0RMAWZm5gBA8XlCiTVIkqQOFJlZzooj5gP/CJyfmasj4kvA88CnMnP6fvM9m5lDrsOJiCXAEoDu7u63Llu2rJQ6O9Xg4CBTp05tdxmVYk8P37Wrd7Dh2X30HDeJq885esj0idjT53a1u4KDm7x74vV0vJuI++l4dsEFF6zNzPljXU+Z1+D0A/2ZuboYvp369TabI2JWZg5ExCxgy3ALZ+ZSYClAT09P9vb2llhq56nVatjT5rKnh+8rG1bBs9uYPn06vb3nDZk+EXs63q/BmTow8Xo63k3E/bQTlHaKKjP/H/CLiOgpRi0AfgosBxYX4xYDd5dVgyRJ6kxl30X1KeDW4g6qnwMfpR6qlkXE5cAm4JKSa5AkSR2m1ICTmQ8Bw51HW1DmdiVJUmfzScaSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyDDiSJKlyjihz5RGxEXgB2Avsycz5EXE88C1gDrAR+GBmPltmHZIkqbO04gjOBZk5LzPnF8NXASsz8zRgZTEsSZLUNO04RbUQ6Cu+9wGL2lCDJEmqsLIDTgL3RsTaiFhSjJuZmQMAxecJJdcgSZI6TKnX4ADnZ+YzEXECsCIiHm90wSIQLQHo7u6mVquVVGJnGhwctKdNZk8P3/btO4rP7cP2biL2dN+udldwcIO7J15Px7uJuJ92glIDTmY+U3xuiYi7gLOBzRExKzMHImIWsGWEZZcCSwF6enqyt7e3zFI7Tq1Ww542lz09fF/ZsAqe3cb06dPp7T1vyPSJ2NN7NrS7goObOjDxejreTcT9tBOUdooqIqZExLSXvwO/DzwKLAcWF7MtBu4uqwZJktSZyjyCMxO4KyJe3s7fZOZ3I+J+YFlEXA5sAi4psQZJktSBSgs4mflz4M3DjN8KLChru5IkST7JWJIkVY4BR5IkVY4BR5IkVY4BR5IkVY4BR5IkVU5DAScizm9knCRJ0njQ6BGcGxocJ0mS1HYHfQ5ORJwHvA3ojoh/t9+k1wKTyyxMkiRptA71oL+jgKnFfNP2G/888IGyipIkSRqLgwaczPwB8IOI+EZmPt2imiS1WateGLl1x28+h9vmvl1Dx7+3p/y6JE18jb6q4TURsRSYs/8ymfmOMoqSJEkai0YDzreBrwI3AXvLK0eSJGnsGg04ezLzK6VWIkmS1CSN3iZ+T0R8IiJmRcTxL/+UWpkkSdIoNXoEZ3HxeeV+4xI4pbnlSJIkjV1DAScz55ZdiCRJUrM0FHAi4rLhxmfmzc0tR5IkaewaPUV11n7fu4AFwAOAAUeSJI07jZ6i+tT+wxFxLHBLKRVJkiSNUaN3UR3oJeC0ZhYiSZLULI1eg3MP9bumoP6SzTcAyxpcdjKwBvjnzLy4uL38W9SfirwR+GBmPnt4ZUvqVK16jYSkia3Ra3Cu3+/7HuDpzOxvcNkrgPXU30AOcBWwMjOvi4iriuHPNrguSZKkQ2roFFXx0s3Hqb9R/Djg140sFxGzgfdQf8XDyxYCfcX3PmBRg7VKkiQ1pNFTVB8E/gKoAQHcEBFXZubth1j0L4H/QD0YvWxmZg4AZOZARJwwwjaXAEsAuru7qdVqjZSqBg0ODtrTJqtST/ftatGGdhavE9+5nX0ba0On7xocfrxGbXB3dfbT8aJKf/tV0ugpqv8InJWZWwAiohv4B2DEgBMRFwNbMnNtRPQebmGZuRRYCtDT05O9vYe9Ch1ErVbDnjZXlXrasutc1q8CtkHXdCbNOW/I5H0ba0ya09uiYjrD1IHq7KfjRZX+9quk0YAz6eVwU9jKoU9vnQ+8LyLeTf3ZOa+NiP8FbI6IWcXRm1nAloOuRZIk6TA1epv4dyPi7yPiIxHxEeDvgO8cbIHMvDozZ2fmHOBDwPcy84+B5fzm3VaLgbtHVbkkSdIIDnoEJyJOpX7NzJUR8X7g7dSvwVkF3DrKbV4HLIuIy4FNwCWjXI8kSdKwDnWK6i+BzwFk5p3AnQARMb+Y9t5GNpKZNeoXKJOZW6m/6kGSJKkUhzpFNScz1x04MjPXUH9QnyRJ0rhzqIDTdZBpRzezEEmSpGY5VMC5PyL+zYEji+tn1pZTkiRJ0tgc6hqcTwN3RcQf8ZtAMx84CvjDEuuSJEkatYMGnMzcDLwtIi4ATi9G/11mfq/0yiRJkkapoQf9Zeb3ge+XXIskSVJTNPqgP0mSpAnDgCNJkirHgCNJkirHgCNJkirHgCNJkirHgCNJkirHgCNJkirHgCNJkirHgCNJkirHgCNJkirHgCNJkirHgCNJkirHgCNJkiqntIATEV0RcV9EPBwRj0XEfy7GHx8RKyLiyeLzuLJqkCRJnanMIzi7gHdk5puBecC7IuJc4CpgZWaeBqwshiVJkpqmtICTdYPF4JHFTwILgb5ifB+wqKwaJElSZzqizJVHxGRgLXAq8OXMXB0RMzNzACAzByLihBGWXQIsAeju7qZWq5VZascZHBy0p01WpZ7u29WiDe3cUXxuZ9/G2tDpuwaHH69RG9xdnf10vKjS336VlBpwMnMvMC8ipgN3RcTph7HsUmApQE9PT/b29pZSY6eq1WrY0+aqUk/v2dCiDa1fBWyDrulMmnPekMn7NtaYNKe3RcV0hqkD1dlPx4sq/e1XSUvuosrM7UANeBewOSJmARSfW1pRgyRJ6hxl3kXVXRy5ISKOBi4EHgeWA4uL2RYDd5dVgyRJ6kxlnqKaBfQV1+FMApZl5t9GxCpgWURcDmwCLimxBkmS1IFKCziZuQ44c5jxW4EFZW1Xmghado2LJHUon2QsSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqx4AjSZIqp7SAExEnR8T3I2J9RDwWEVcU44+PiBUR8WTxeVxZNUiSpM5U5hGcPcC/z8w3AOcC/zYi3ghcBazMzNOAlcWwJElS05QWcDJzIDMfKL6/AKwHTgIWAn3FbH3AorJqkCRJnemIVmwkIuYAZwKrgZmZOQD1EBQRJ4ywzBJgCUB3dze1Wq0VpXaMwcFBe9pkh9PTfbvKrWXC2Lmj+NzOvo21odN3DQ4/XqM2uNu//Wbz39PxqfSAExFTgTuAT2fm8xHR0HKZuRRYCtDT05O9vb2l1diJarUa9rS5Dqen92wot5YJY/0qYBt0TWfSnPOGTN63scakOb0tL6vKpg74t99s/ns6PpV6F1VEHEk93NyamXcWozdHxKxi+ixgS5k1SJKkzlPmXVQB/BWwPjP/236TlgOLi++LgbvLqkGSJHWmMk9RnQ98GHgkIh4qxn0OuA5YFhGXA5uAS0qsQZIkdaDSAk5m/hgY6YKbBWVtV5IkyScZS5KkyjHgSJKkymnJc3Ckdmj1rdj7dnn7tySNFx7BkSRJlWPAkSRJlWPAkSRJlWPAkSRJlWPAkSRJlWPAkSRJlWPAkSRJlWPAkSRJlWPAkSRJlWPAkSRJlWPAkSRJlWPAkSRJlWPAkSRJlWPAkSRJlWPAkSRJlWPAkSRJlVNawImIr0fEloh4dL9xx0fEioh4svg8rqztS5KkzlXmEZxvAO86YNxVwMrMPA1YWQxLkiQ1VWkBJzN/CGw7YPRCoK/43gcsKmv7kiSpcx3R4u3NzMwBgMwciIgTRpoxIpYASwC6u7up1WqtqbBDDA4OVr6n+3a1eIO7Btm3sdbijU5wO3cUn9uH7509bbrB3dX/22+1Tvj3dCJqdcBpWGYuBZYC9PT0ZG9vb3sLqpharUbVe3rPhtZub9/GGpPm9LZ2oxPd+lXANuiazqQ55w2ZbE+bb+pA9f/2W60T/j2diFp9F9XmiJgFUHxuafH2JUlSB2h1wFkOLC6+LwbubvH2JUlSByjzNvHbgFVAT0T0R8TlwHXAOyPiSeCdxbAkSVJTlXYNTmZeOsKkBWVtU5IkCXySsSRJqiADjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqhwDjiRJqpxx+7JNjX+tfpmlJEmN8giOJEmqHAOOJEmqHAOOJEmqHAOOJEmqHC8yHsfKvIh33y4vEpYkVZdHcCRJUuUYcCRJUuUYcCRJUuUYcCRJUuUYcCRJUuW0JeBExLsiYkNE/CwirmpHDZIkqbpaHnAiYjLwZeAPgDcCl0bEG1tdhyRJqq52HME5G/hZZv48M38NfBNY2IY6JElSRbXjQX8nAb/Yb7gfOOfAmSJiCbCkGNwVEY+2oLZO8jrgV+0uomLs6Sg9Daz8s2En2dPms6fNZ0+bq6cZK2lHwIlhxuWQEZlLgaUAEbEmM+eXXVgnsafNZ0+bz542nz1tPnvaXBGxphnraccpqn7g5P2GZwPPtKEOSZJUUe0IOPcDp0XE3Ig4CvgQsLwNdUiSpIpq+SmqzNwTEZ8E/h6YDHw9Mx87xGJLy6+s49jT5rOnzWdPm8+eNp89ba6m9DMyh1z+IkmSNKH5JGNJklQ5BhxJklQ5bQ04h3plQ0QcGxH3RMTDEfFYRHy00WU71Rh7ujEiHomIh5p1m14VNNDT4yLirohYFxH3RcTpjS7bqcbYU/fTYUTE1yNiy0jPDIu6/170fF1EvGW/ae6nBxhjP91Hh9FAT18fEasiYldEfOaAaYe/j2ZmW36oX2D8T8ApwFHAw8AbD5jnc8CfF9+7gW3FvIdcthN/xtLTYngj8Lp2/x7j6afBnv4F8Pni++uBlY0u24k/Y+lpMex+Onxffw94C/DoCNPfDfwf6s8iOxdY3ej/Hp34M9p+FtPcR0fX0xOAs4D/Cnxmv/Gj2kfbeQSnkVc2JDAtIgKYSv0/xnsaXLYTjaWnGl4jPX0jsBIgMx8H5kTEzAaX7URj6alGkJk/pP73PJKFwM1Z94/A9IiYhfvpsMbQT43gUD3NzC2ZeT+w+4BJo9pH2xlwhntlw0kHzHMj8AbqDwJ8BLgiM/c1uGwnGktPoR5+7o2ItcWrMtRYTx8G3g8QEWcDv039AZbup8MbS0/B/XS0Ruq7++noHKxv7qPNNap9tB2vanhZI69suAh4CHgH8DvAioj4UYPLdqJR9zQznwfOz8xnIuKEYvzjReLuZI309DrgSxHxEPXQ+CD1o2Lup8MbS0/B/XS0Ruq7++noHKxv7qPNNap9tJ1HcBp5ZcNHgTuLQ4A/A56ifj7e1z0Mbyw9JTOfKT63AHdRPyzY6Q7Z08x8PjM/mpnzgMuoX9v0VCPLdqix9NT9dPRG6rv76eiM2Df30aYb1T7azoDTyCsbNgELAIrz7z3AzxtcthONuqcRMSUiphXjpwC/D/gG9wZ6GhHTi2kA/xr4YXFEzP10eKPuqfvpmCwHLivu/jkXeC4zB3A/Ha1h++k+WopR7aNtO0WVI7yyISL+pJj+VeC/AN+IiEeoH6L6bGb+CmC4Zdvxe4wnY+lpRJwC3FW/9pgjgL/JzO+25RcZRxrs6RuAmyNiL/BT4PKDLduO32M8GUtPgZm4nw4rIm4DeoHXRUQ/8HngSHilp9+hfufPz4CXqB/NdT8dwWj7ifvoiA7V04j4LWAN8FpgX0R8mvrdUs+PZh/1VQ2SJKlyfJKxJEmqHAOOJEmqHAOOJEmqHAOOJEmqHAOOJEmqHAOOpKaIiD+MiIyI17e7Fkky4EhqlkuBH1N/CFcpImJyWeuWVC0GHEljFhFTgfOpP5DvQ8W4yRFxfUQ8EhHrIuJTxfizIuInEfFwRNwXEdMi4iMRceN+6/vbiOgtvg9GxBcjYjVwXkT8p4i4PyIejYilUTxRLSJOjYh/KNb7QET8TkTcEhEL91vvrRHxvlb1RVL7GHAkNcMi4LuZ+QSwLSLeAiwB5gJnZuYZwK3FY9a/Rf0t9m8GLgR2HGLdU4BHM/OczPwxcGNmnpWZpwNHAxcX890KfLlY79uAAeAmiifMRsSxxfjvNOuXljR+GXAkNcOlwDeL798shi8EvpqZewAycxv1d58NZOb9xbjnX55+EHuBO/YbviAiVhevG3kH8LvFu39Oysy7ivXuzMyXMvMHwKnFW50vBe5oYHuSKqBt76KSVA0RMYN60Dg9IpL6u2ISWFt8vmr2YcYB7OHV/4era7/vOzNzb7GtLuB/APMz8xcR8YVi3jhIibcAf0T91NnHGvy1JE1wHsGRNFYfAG7OzN/OzDmZeTLwFPAA8CcRcQRARBwPPA6cGBFnFeOmFdM3AvMiYlJEnAycPcK2Xg4+vyqu+/kA1I8EAf0RsahY72si4phi3m8Any7m6/iXSEqdwoAjaawuBe46YNwdwInAJmBdRDwM/KvM/DXwL4EbinErqIeW/0s9FD0CXE89HA2RmduB/1nM97+B+/eb/GHgTyNiHfAT4LeKZTYD64G/HuPvKWkC8W3ikiqtOJLzCPCWzHyu3fVIag2P4EiqrIi4kPppsRsMN1Jn8QiOJEmqHI/gSJKkyjHgSJKkyjHgSJKkyjHgSJKkyjHgSJKkyvn/g07z/6HVZnkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "\n",
    "ax.vlines(bootstrap_train_mean, [0], 80, lw=2.5, linestyle=\"-\", label=\"Mean\")\n",
    "\n",
    "ax.hist(\n",
    "    bootstrap_train_accuracies, bins=7, color=\"#0080ff\", edgecolor=\"none\", alpha=0.3\n",
    ")\n",
    "plt.legend(loc=\"upper left\")\n",
    "\n",
    "plt.xlim([0.8, 1.1])\n",
    "\n",
    "plt.xlabel('Accuracy')\n",
    "plt.ylabel('Count')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.grid()\n",
    "plt.savefig(\"matplotlib-figures/bootstrap-hist.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Useful tidbit: Here, we are creating validation sets from the training set. However, if you don't tune your model on the training set, you could use the whole dataset for this since actually don't need a test set for evaluation then. Or in other words, the bootstrap methods can serve can take the place of a test set for model evaluation. This is particularly attractive for small datasets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1) A *t* Confidence Interval Interval from Bootstrap Samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Assuming that the sample means are normal distributed we could compute the confidence interval formula as before, as follows:\n",
    "\n",
    "$$\n",
    "\\text{ACC}_{\\text{test}} \\pm z \\times \\text{SE}.\n",
    "$$\n",
    "\n",
    "- But since we are dealing with a relatively small sample where we also estimate the population standard deviation via the sample standard deviation (this is used t calculate the standard error), we replace the $z$ value with a $t$ value from a t-distribution:\n",
    "\n",
    "$$\n",
    "\\text{ACC}_{\\text{test}} \\pm t \\times \\text{SE}.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- We can then compute the standard error (SE) by the as the standard deviation of this distribution:\n",
    "\n",
    "$$\n",
    "\\text{SE}=\\sqrt{\\frac{1}{b-1} \\sum_{j=1}^{b}\\left(\\text{ACC}_{\\text{boot},j}-{\\text{ACC}_\\text{bootavg}}\\right)^{2}}.\n",
    "$$\n",
    "\n",
    "\n",
    "- Note that we usually divide the standard deviation (SD) by $\\sqrt{n}$ to obtain the standard error (SE) where $n = b$:\n",
    "\n",
    "$$ \\text{SE} = \\frac{\\text{SD}}{\\sqrt{n}}.$$\n",
    "\n",
    "- However, this is not necessary here since the bootstrap distribution is a distribution of means (as opposed to single data points), and $\\text{ACC}_{\\text{bootavg}}$ is the mean of these means.\n",
    "\n",
    "- (As an optional exercise, you can try to modify the code below to include the division by $\\sqrt{n}$ (where $\\sqrt{n} = \\sqrt{b}$), and you will probably find that this shrinks the confidence interval to an unrealistic degree, which also doesn't match the percentile method results anymore that we will introduce shortly.)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Again, instead of using a t-table in our old statistics textbook, let's use SciPy to obtain the t-value for a 95% confidence interval and $b-1$ degrees of freedom:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.971956544249395\n"
     ]
    }
   ],
   "source": [
    "confidence = 0.95  # Change to your desired confidence level\n",
    "t_value = scipy.stats.t.ppf((1 + confidence) / 2.0, df=bootstrap_rounds - 1)\n",
    "print(t_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.879470830164037 1.0132047668825668\n"
     ]
    }
   ],
   "source": [
    "se = 0.0\n",
    "for acc in bootstrap_train_accuracies:\n",
    "    se += (acc - bootstrap_train_mean) ** 2\n",
    "se = np.sqrt((1.0 / (bootstrap_rounds - 1)) * se)\n",
    "\n",
    "ci_length = t_value * se\n",
    "\n",
    "ci_lower = bootstrap_train_mean - ci_length\n",
    "ci_upper = bootstrap_train_mean + ci_length\n",
    "\n",
    "print(ci_lower, ci_upper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Now, let's add the CI to our histogram plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEYCAYAAABRMYxdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhoklEQVR4nO3de5hddXno8e+bG5ObJIEhTYhtQpNnEEMMGgKIlkG0UFFJVVRO0TRG09aDd5CL5YiX80AVz4OChxpodeRYWgpEQosUGhmgGMGAEILhKjEGpolJCGTIxYS854/ZxCQzSXZm9p49s/b38zzz7L1+v3V59+ta4XXdfpGZSJIkFcmAWgcgSZJUaRY4kiSpcCxwJElS4VjgSJKkwrHAkSRJhTOo1gGUY9SoUTl58uRah1EoL7/8MsOHD691GIViTivPnFaeOa08c1pZDz744NrMbOzpevpFgTN27FiWLFlS6zAKpbW1lebm5lqHUSjmtPLMaeWZ08ozp5UVEb+uxHq8RCVJkgrHAkeSJBWOBY4kSSqcfnEPTle2bdvGqlWr2LJlS61D6VMaGhqYMGECgwcPrnUokiTVTL8tcFatWsXIkSOZOHEiEVHrcPqEzGTdunWsWrWKSZMm1TocSZJqpt9eotqyZQuHHHKIxc0uIoJDDjnEs1qSpLpX1QInIj4bEY9FxLKIuD4iGiJiTETcGRFPlT5H92D9lQy3EMyJJElVLHAi4nDgU8CMzJwKDAQ+BFwALMrMKcCi0rQkSVLFVPsS1SBgaEQMAoYBzwNnAC2l/hZgVpVjqJqI4MMf/vDO6e3bt9PY2Mi73vWuGkYlSZKqVuBk5nPA5cBKoA14MTPvAMZmZltpnjbgsGrFUG3Dhw9n2bJlbN68GYA777yTww8/vMZRSZKkqj1FVbq35gxgErAB+NeIOPsAlp8HzANobGyktbV1t/6DDz6YjRs3VircbjvllFO48cYbmTVrFj/4wQ9473vfy09/+lM2btzIyy+/zHnnncdjjz3GK6+8woUXXsjpp5/Or3/9a+bNm8emTZsAuPzyyznuuOO49957ufTSSznkkEP45S9/yfTp07n22msP+L6aLVu2dMrXntrb2/c7jw6MOT1w967axtrNyaFDg7dO6PxqA3Naeea08sxp31TNx8TfDjybmb8FiIibgTcDqyNiXGa2RcQ4YE1XC2fmfGA+QFNTU+45zsfy5csZOXLkzul/XfIbbnxwVcWCf/+bJnDmjNfud76PfOQjfOUrX+HMM89k+fLl/NVf/RUPPPAAI0eO5NJLL+XUU0/luuuuY8OGDcycOZN3v/vdHHHEEfzkJz+hoaGBp556irPOOoslS5YwbNgwli5dymOPPcb48eM58cQTWbp0KW95y1sOKPaGhgaOOeaYfc7j2CmVZ04P3NXfXcz9z67nuEljuLj5hE795rTyzGnlmdO+qZoFzkrg+IgYBmwGTgGWAC8Ds4HLSp+3VGJjq17YzP3Prq/EqgA4/ohDyppv2rRprFixguuvv553vvOdu/XdcccdLFy4kMsvvxzoOLOycuVKxo8fzznnnMPDDz/MwIEDefLJJ3cuM3PmTCZMmADA9OnTWbFixQEXOJIk1buqFTiZeX9E3Ag8BGwHfkHHGZkRwA0RMZeOIujMSmxvwuihHDdpTCVWtXN95XrPe97DueeeS2trK+vWrdvZnpncdNNNNDU17Tb/JZdcwtixY3nkkUfYsWMHDQ0NO/sOOuignd8HDhzI9u3be/ArJEmqT1V9k3Fmfgn40h7NW+k4m1NRZ854bVmXlKrhox/9KAcffDBHH330btdhTz31VK688kquvPJKIoJf/OIXHHPMMbz44otMmDCBAQMG0NLSwiuvvFKTuCVJKqp++ybjvmTChAl8+tOf7tR+8cUXs23bNqZNm8bUqVO5+OKLAfjEJz5BS0sLxx9/PE8++STDhw/v7ZAlSSq0fjsWVV/Q3t7eqa25uXnnzWZDhw7lu9/9bqd5pkyZwtKlS3dOX3rppZ2WBbjqqqsqG7AkSXXCMziSJKlwLHAkSVLhWOBIkqTCscCRJEmFY4EjSZIKxwJHkiQVjgVOD3zrW99i6tSpvP71r+eKK67Y2X7JJZdw+OGHM336dKZPn85tt90GwH333ce0adM49thjefrppwHYsGEDp556KpnZ5Ta2bdvGBRdcwJQpU5g6dSozZ87kxz/+MQATJ05k7dq11f2RkiT1Q74Hp5uWLVvGNddcwwMPPMCQIUM47bTTOP3005kyZQoAn/3sZzn33HN3W+ab3/wmN910EytWrODqq6/mm9/8Jl/96le56KKL9jpi+MUXX0xbWxvLli3joIMOYvXq1dx9991V/32SJPVnnsHppuXLl3P88cczbNgwBg0axEknncSCBQv2uczgwYPZvHkzmzZtYvDgwTzzzDM899xznHTSSV3Ov2nTJq655hquvPLKnWNUjR07lg984AMV/z2SJBVJYQqcHz39I+bcPoc5t8/p1Hfe3ecx5/Y5XPvotbu1P77+8Z3LPL7+8d361m7e96WfqVOncs8997Bu3To2bdrEbbfdxm9+85ud/VdddRXTpk3jox/9KC+88AIAF154IfPmzeOKK67gnHPO4Ytf/CJf/epX97qNp59+mj/8wz/kNa95zX5/vyRJ+r3CFDjPtz/PktVLWLJ6Sae+pb9dypLVS3j2xWd3a9/4u407l9n4u4279W19Zes+t/e6172O888/n3e84x2cdtppvOENb2DQoI4rfn/zN3/DM888w8MPP8y4ceP4/Oc/D8D06dP52c9+xl133cWvfvUrxo8fT2bywQ9+kLPPPpvVq1f3JAWSJKmkMAXO+BHjmTF2BjPGzujUN61xGjPGzmDSwZN2ax85ZOTOZUYOGblb30EDD9rvNufOnctDDz3EPffcw5gxY3befzN27FgGDhzIgAED+PjHP84DDzyw23KZyde+9jUuvvhivvzlL/PlL3+Zs88+m29/+9u7zTd58mRWrlzJxo27F1+SJGnfCnOT8azJs5g1eVaXfd846Rtdth855ki+d9r3uuw7dOih+93mmjVrOOyww1i5ciU333wzixcvBqCtrY1x48YBsGDBAqZOnbrbci0tLZx++umMHj2aTZs2MWDAAAYMGMCmTZt2m2/YsGHMnTuXT33qU3z3u99lyJAhtLW1sWjRIs4+++z9xidJUr0qTIFTC+973/tYt24dgwcP5jvf+Q6jR48G4Atf+AIPP/wwEcHEiRN3G1F806ZNtLS0cMcddwDwuc99jve9730MGTKE66+/vtM2vva1r/G3f/u3HHXUUTQ0NDB8+HC+8pWv9M4PlCSpn7LA6YF77723y/brrrtur8sMGzaMu+66a+f0W9/6Vh599NG9zj9kyBC+/vWv8/Wvf71T34oVK8oPVpKkOlKYe3AkSZJeZYEjSZIKp2oFTkQ0RcTDu/y9FBGfiYgxEXFnRDxV+hzd3W3sbXiDemZOJEmqYoGTmU9k5vTMnA68CdgELAAuABZl5hRgUWn6gDU0NLBu3Tr/g76LzGTdunU0NDTUOhRJkmqqt24yPgV4JjN/HRFnAM2l9hagFTj/QFc4YcIEVq1axW9/+9uKBVkEDQ0NTJgwodZhSJJUU71V4HwIePUZ6LGZ2QaQmW0RcVh3Vjh48GAmTZq0/xklSVLdiWpf4omIIcDzwOszc3VEbMjMUbv0v5CZne7DiYh5wDyAxsbGN91www1VjbPetLe3M2LEiFqHUSjm9MBdev9mnnhhB02jB3DhcUM79ZvTyjOnlWdOK+vkk09+MDM7D0twgHrjDM6fAQ9l5qsDLa2OiHGlszfjgDVdLZSZ84H5AE1NTdnc3NwLodaP1tZWzGllmdMDd/UTi+GF9YwaNYrm5hM69ZvTyjOnlWdO+6beeEz8LH5/eQpgITC79H02cEsvxCBJkupIVQuciBgGvAO4eZfmy4B3RMRTpb7LqhmDJEmqP1W9RJWZm4BD9mhbR8dTVZIkSVXhm4wlSVLhWOBIkqTCscCRJEmFY4EjSZIKxwJHkiQVjgWOJEkqHAscSZJUOBY4kiSpcCxwJElS4VjgSJKkwrHAkSRJhWOBI0mSCscCR5IkFY4FjiRJKhwLHEmSVDgWOJIkqXAscCRJUuFY4EiSpMKxwJEkSYVT1QInIkZFxI0R8XhELI+IEyJiTETcGRFPlT5HVzMGSZJUfwZVef3fAm7PzPdHxBBgGHARsCgzL4uIC4ALgPOrHIck9Ypbn6h1BPs2stYBSL2kamdwIuI1wJ8A/wCQmb/LzA3AGUBLabYWYFa1YpAkSfWpmpeojgB+C3wvIn4REddGxHBgbGa2AZQ+D6tiDJIkqQ5FZlZnxREzgJ8BJ2bm/RHxLeAl4JOZOWqX+V7IzE734UTEPGAeQGNj45tuuOGGqsRZr9rb2xkxYkStwygUc3rgLr1/M0+8sIOm0QO48Lihnfr7Y05f3FrrCPZt4Lb+l9O+rj/up33ZySef/GBmzujpeqp5D84qYFVm3l+avpGO+21WR8S4zGyLiHHAmq4Wzsz5wHyApqambG5urmKo9ae1tRVzWlnm9MBd/cRieGE9o0aNorn5hE79/TGnff0enBFt/S+nfV1/3E/rQdUuUWXmfwO/iYimUtMpwC+BhcDsUtts4JZqxSBJkupTtZ+i+iTww9ITVL8C5tBRVN0QEXOBlcCZVY5BkiTVmaoWOJn5MNDVdbRTqrldSZJU33yTsSRJKhwLHEmSVDgWOJIkqXAscCRJUuFY4EiSpMKxwJEkSYVjgSNJkgrHAkeSJBWOBY4kSSocCxxJklQ4FjiSJKlwLHAkSVLhWOBIkqTCscCRJEmFY4EjSZIKxwJHkiQVjgWOJEkqHAscSZJUOBY4kiSpcAZVc+URsQLYCLwCbM/MGRExBvgXYCKwAvhAZr5QzTgkSVJ96Y0zOCdn5vTMnFGavgBYlJlTgEWlaUmSpIqpxSWqM4CW0vcWYFYNYpAkSQVW7QIngTsi4sGImFdqG5uZbQClz8OqHIMkSaozVb0HBzgxM5+PiMOAOyPi8XIXLBVE8wAaGxtpbW2tUoj1qb293ZxWmDk9cBs2bC59bugyd/0xpzu21jqCfWvf1v9y2tf1x/20HlS1wMnM50ufayJiATATWB0R4zKzLSLGAWv2sux8YD5AU1NTNjc3VzPUutPa2oo5rSxzeuCufmIxvLCeUaNG0dx8Qqf+/pjTW5+odQT7NqKt/+W0r+uP+2k9qNolqogYHhEjX/0O/CmwDFgIzC7NNhu4pVoxSJKk+lTNMzhjgQUR8ep2/ikzb4+InwM3RMRcYCVwZhVjkCRJdahqBU5m/gp4Qxft64BTqrVdSZIk32QsSZIKxwJHkiQVjgWOJEkqHAscSZJUOBY4kiSpcMoqcCLixHLaJEmS+oJyz+BcWWabJElSze3zPTgRcQLwZqAxIj63S9drgIHVDEySJKm79veivyHAiNJ8I3dpfwl4f7WCkiRJ6ol9FjiZeTdwd0R8PzN/3UsxSaqx3howct3m3392tc0dWzu3v7up+nFJ6v/KHarhoIiYD0zcdZnMfFs1gpIkSeqJcgucfwX+HrgWeKV64UiSJPVcuQXO9sy8uqqRSJIkVUi5j4nfGhGfiIhxETHm1b+qRiZJktRN5Z7BmV36PG+XtgSOqGw4kiRJPVdWgZOZk6odiCRJUqWUVeBExEe6as/MH1Q2HEmSpJ4r9xLVsbt8bwBOAR4CLHAkSVKfU+4lqk/uOh0RBwPXVSUiSZKkHir3Kao9bQKmVDIQSZKkSin3Hpxb6XhqCjoG2XwdcEOZyw4ElgDPZea7So+X/wsdb0VeAXwgM184sLAl1aveGkZCUv9W7j04l+/yfTvw68xcVeaynwaW0zECOcAFwKLMvCwiLihNn1/muiRJkvarrEtUpUE3H6djRPHRwO/KWS4iJgCn0zHEw6vOAFpK31uAWWXGKkmSVJZyL1F9APgG0AoEcGVEnJeZN+5n0SuAL9BRGL1qbGa2AWRmW0QctpdtzgPmATQ2NtLa2lpOqCpTe3u7Oa2wIuV0x9Ze2tCW0nDiWzawY0Vr5/6t7V23q9vatxVnP+0rinTsF0m5l6i+CBybmWsAIqIR+E9grwVORLwLWJOZD0ZE84EGlpnzgfkATU1N2dx8wKvQPrS2tmJOK6tIOe21+1yWLwbWQ8MoBkw8oVP3jhWtDJjY3EvB1IcRbcXZT/uKIh37RVJugTPg1eKmZB37v7x1IvCeiHgnHe/OeU1E/D9gdUSMK529GQes2edaJEmSDlC5j4nfHhH/ERF/GRF/Cfw7cNu+FsjMCzNzQmZOBD4E/CQzzwYW8vuxrWYDt3QrckmSpL3Y5xmciJhMxz0z50XEe4G30HEPzmLgh93c5mXADRExF1gJnNnN9UiSJHVpf5eorgAuAsjMm4GbASJiRqnv3eVsJDNb6bhBmcxcR8dQD5IkSVWxv0tUEzNz6Z6NmbmEjhf1SZIk9Tn7K3Aa9tE3tJKBSJIkVcr+CpyfR8TH92ws3T/zYHVCkiRJ6pn93YPzGWBBRPwFvy9oZgBDgD+vYlySJEndts8CJzNXA2+OiJOBqaXmf8/Mn1Q9MkmSpG4q60V/mXkXcFeVY5EkSaqIcl/0J0mS1G9Y4EiSpMKxwJEkSYVjgSNJkgrHAkeSJBWOBY4kSSocCxxJklQ4FjiSJKlwLHAkSVLhWOBIkqTCscCRJEmFY4EjSZIKxwJHkiQVTtUKnIhoiIgHIuKRiHgsIr5cah8TEXdGxFOlz9HVikGSJNWnap7B2Qq8LTPfAEwHTouI44ELgEWZOQVYVJqWJEmqmKoVONmhvTQ5uPSXwBlAS6m9BZhVrRgkSVJ9GlTNlUfEQOBBYDLwncy8PyLGZmYbQGa2RcRhe1l2HjAPoLGxkdbW1mqGWnfa29vNaYUVKac7tvbShrZsLn1uYMeK1s79W9u7ble3tW8rzn7aVxTp2C+SqhY4mfkKMD0iRgELImLqASw7H5gP0NTUlM3NzVWJsV61trZiTiurSDm99Yle2tDyxcB6aBjFgIkndOresaKVARObeymY+jCirTj7aV9RpGO/SHrlKarM3AC0AqcBqyNiHEDpc01vxCBJkupHNZ+iaiyduSEihgJvBx4HFgKzS7PNBm6pVgySJKk+VfMS1TigpXQfzgDghsz8t4hYDNwQEXOBlcCZVYxBkiTVoaoVOJm5FDimi/Z1wCnV2q7UH/TaPS6SVKd8k7EkSSocCxxJklQ4FjiSJKlwLHAkSVLhWOBIkqTCscCRJEmFY4EjSZIKxwJHkiQVjgWOJEkqHAscSZJUOBY4kiSpcCxwJElS4VjgSJKkwrHAkSRJhWOBI0mSCscCR5IkFY4FjiRJKhwLHEmSVDgWOJIkqXCqVuBExGsj4q6IWB4Rj0XEp0vtYyLizoh4qvQ5uloxSJKk+lTNMzjbgc9n5uuA44H/GRFHARcAizJzCrCoNC1JklQxVStwMrMtMx8qfd8ILAcOB84AWkqztQCzqhWDJEmqT4N6YyMRMRE4BrgfGJuZbdBRBEXEYXtZZh4wD6CxsZHW1tbeCLVutLe3m9MKO5Cc7tha3Vj6jS2bS58b2LGitXP/1vau29Vt7ds89ivNf0/7pqoXOBExArgJ+ExmvhQRZS2XmfOB+QBNTU3Z3NxctRjrUWtrK+a0sg4kp7c+Ud1Y+o3li4H10DCKARNP6NS9Y0UrAyY293pYRTaizWO/0vz3tG+q6lNUETGYjuLmh5l5c6l5dUSMK/WPA9ZUMwZJklR/qvkUVQD/ACzPzP+zS9dCYHbp+2zglmrFIEmS6lM1L1GdCHwYeDQiHi61XQRcBtwQEXOBlcCZVYxBkiTVoaoVOJn5X8Debrg5pVrblSRJ8k3GkiSpcCxwJElS4fTKe3CkWujtR7F3bPXxb0nqKzyDI0mSCscCR5IkFY4FjiRJKhwLHEmSVDgWOJIkqXB8ikr91o+e/hHPtz/P+BHjmTV5Vq3DkeqWx6L6Igsc9Vu3PH0LS1YvYcbYGf6jKtWQx6L6Igsc9VtHjjlyt09JteGxqL7IAkf91vkzz691CJLwWFTf5E3GkiSpcCxwJElS4XiJSv3W2s1r2frKVg4aeBCHDj201uFIdctjUX2RBY76rfPuPm/nkxvfO+17tQ5Hqlsei+qLvEQlSZIKxzM46rfmHj2XMyaf4SlxqcY8FtUXWeCo33rL4W+pdQiS8FhU31S1S1QR8Y8RsSYilu3SNiYi7oyIp0qfo6u1fUmSVL+qeQ/O94HT9mi7AFiUmVOARaVpSZKkiqraJarMvCciJu7RfAbQXPreArQCvgJT3eIAf1Lf4LGovqi378EZm5ltAJnZFhGH7W3GiJgHzANobGyktbW1dyKsE+3t7f0+py3/3cLTW59m8kGTGbVqVKf+HVt7OaCt7exY0drLG+3ntmwufW7oOnfmtOLat1X+2N/fsVh0Rfj3tIj67E3GmTkfmA/Q1NSUzc3NtQ2oYFpbW+nvOW25vQVWw6hRo7r8Lbc+0bvx7FjRyoCJnePQPixfDKyHhlEMmHhCp25zWnkj2ip/7O/vWCy6Ivx7WkS9XeCsjohxpbM344A1vbx9FYgvFJP6Bo9F9UW9/aK/hcDs0vfZwC29vH1JklQHqvmY+PXAYqApIlZFxFzgMuAdEfEU8I7StCRJUkVV8ymqs/bSdUq1tilJkgSORaV+7Ly7z+PUG0/lvLvPq3UoUl3zWFRfZIGjfmvt5rU8//LzrN28ttahSHXNY1F9UZ99TFzanxMPP5HxI8Yz6eBJtQ5Fqmsei+qLLHDUb33s6I/VOgRJeCyqb/ISlSRJKhwLHEmSVDheolK/9fj6x9n4u42MHDKSI8ccWetwpLrlsai+yAJH/dbfPfB3LFm9hBljZ/iqeKmGPBbVF1ngqNt6ezDLPa3b9PvPWsciSepbLHDUb50+8Xy2bN9Iw6CRtQ5Fqmvnzzx/5yUqqa+wwFG/NX641/qlvsD7btQX+RSVJEkqHAscSZJUOF6i6sOqeePsjq39/8bcu5+7lt9ufpbGoZM46XDfpCrVyrWPXsuzLz7LpIMn+VZj9RmewVG/9eSG+3ho7UKe3HBfrUOR6tp9z93HwmcWct9zHovqOyxw1G+NGHIoo4aMZ8SQQ2sdilTXDh16KOOHj+fQoR6L6ju8RKV+66wp36h1CJKAb5zksai+xzM4kiSpcCxwJElS4dSkwImI0yLiiYh4OiIuqEUMkiSpuHq9wImIgcB3gD8DjgLOioijejsO9X/XPDaHi352NNc8NqfWoUh1bc7tczi65Wjm3O6xqL6jFmdwZgJPZ+avMvN3wD8DZ9QgDkmSVFCRmb27wYj3A6dl5sdK0x8GjsvMc/aYbx4wrzQ5FVjWq4EW36HA2loHUTDmtPLMaeWZ08ozp5XVlJk9Hrm1Fo+JRxdtnaqszJwPzAeIiCWZOaPagdUTc1p55rTyzGnlmdPKM6eVFRFLKrGeWlyiWgW8dpfpCcDzNYhDkiQVVC0KnJ8DUyJiUkQMAT4ELKxBHJIkqaB6/RJVZm6PiHOA/wAGAv+YmY/tZ7H51Y+s7pjTyjOnlWdOK8+cVp45rayK5LPXbzKWJEmqNt9kLEmSCscCR5IkFU5NC5z9DdkQEQdHxK0R8UhEPBYRc8pdtl71MKcrIuLRiHi4Uo/pFUEZOR0dEQsiYmlEPBARU8tdtl71MKfup12IiH+MiDUR0eU7w6LDt0s5XxoRb9ylz/10Dz3Mp/toF8rI6ZERsTgitkbEuXv0Hfg+mpk1+aPjBuNngCOAIcAjwFF7zHMR8Hel743A+tK8+122Hv96ktPS9Arg0Fr/jr70V2ZOvwF8qfT9SGBRucvW419Pclqadj/tOq9/ArwRWLaX/ncCP6bjXWTHA/eX+79HPf51N5+lPvfR7uX0MOBY4H8D5+7S3q19tJZncMoZsiGBkRERwAg6/mO8vcxl61FPcqqulZPTo4BFAJn5ODAxIsaWuWw96klOtReZeQ8dx/PenAH8IDv8DBgVEeNwP+1SD/KpvdhfTjNzTWb+HNi2R1e39tFaFjiHA7/ZZXpVqW1XVwGvo+NFgI8Cn87MHWUuW496klPoKH7uiIgHS0NlqLycPgK8FyAiZgJ/RMcLLN1Pu9aTnIL7aXftLe/up92zr7y5j1ZWt/bRWgzV8Kpyhmw4FXgYeBvwx8CdEXFvmcvWo27nNDNfAk7MzOcj4rBS++OliruelZPTy4BvRcTDdBSNv6DjrJj7add6klNwP+2uveXd/bR79pU399HK6tY+WsszOOUM2TAHuLl0CvBp4Fk6rsc73EPXepJTMvP50ucaYAEdpwXr3X5zmpkvZeaczJwOfISOe5ueLWfZOtWTnLqfdt/e8u5+2j17zZv7aMV1ax+tZYFTzpANK4FTAErX35uAX5W5bD3qdk4jYnhEjCy1Dwf+FEdwhzJyGhGjSn0AHwPuKZ0Rcz/tWrdz6n7aIwuBj5Se/jkeeDEz23A/7a4u8+k+WhXd2kdrdokq9zJkQ0T8dan/74GvAt+PiEfpOEV1fmauBehq2Vr8jr6kJzmNiCOABR33HjMI+KfMvL0mP6QPKTOnrwN+EBGvAL8E5u5r2Vr8jr6kJzkFxuJ+2qWIuB5oBg6NiFXAl4DBsDOnt9Hx5M/TwCY6zua6n+5Fd/OJ++he7S+nEfEHwBLgNcCOiPgMHU9LvdSdfdShGiRJUuH4JmNJklQ4FjiSJKlwLHAkSVLhWOBIkqTCscCRJEmFY4EjqSIi4s8jIiPiyFrHIkkWOJIq5Szgv+h4CVdVRMTAaq1bUrFY4EjqsYgYAZxIxwv5PlRqGxgRl0fEoxGxNCI+WWo/NiJ+GhGPRMQDETEyIv4yIq7aZX3/FhHNpe/tEfGViLgfOCEi/ldE/DwilkXE/Ci9US0iJkfEf5bW+1BE/HFEXBcRZ+yy3h9GxHt6Ky+SascCR1IlzAJuz8wngfUR8UZgHjAJOCYzpwE/LL1m/V/oGMX+DcDbgc37WfdwYFlmHpeZ/wVclZnHZuZUYCjwrtJ8PwS+U1rvm4E24FpKb5iNiINL7bdV6kdL6rsscCRVwlnAP5e+/3Np+u3A32fmdoDMXE/H2GdtmfnzUttLr/bvwyvATbtMnxwR95eGG3kb8PrS2D+HZ+aC0nq3ZOamzLwbmFwa1fks4KYytiepAGo2FpWkYoiIQ+goNKZGRNIxVkwCD5Y+d5u9izaA7ez+f7gadvm+JTNfKW2rAfi/wIzM/E1EXFKaN/YR4nXAX9Bx6eyjZf4sSf2cZ3Ak9dT7gR9k5h9l5sTMfC3wLPAQ8NcRMQggIsYAjwPjI+LYUtvIUv8KYHpEDIiI1wIz97KtVwuftaX7ft4PHWeCgFURMau03oMiYlhp3u8DnynNV/eDSEr1wgJHUk+dBSzYo+0mYDywElgaEY8A/yMzfwd8ELiy1HYnHUXLfXQURY8Cl9NRHHWSmRuAa0rz/Qj4+S7dHwY+FRFLgZ8Cf1BaZjWwHPheD3+npH7E0cQlFVrpTM6jwBsz88VaxyOpd3gGR1JhRcTb6bgsdqXFjVRfPIMjSZIKxzM4kiSpcCxwJElS4VjgSJKkwrHAkSRJhWOBI0mSCuf/A5wgChz+J6JKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "ax.vlines(bootstrap_train_mean, [0], 80, lw=2.5, linestyle=\"-\", label=\"Mean\")\n",
    "\n",
    "ax.vlines(ci_lower, [0], 15, lw=2.5, linestyle=\"dotted\", label=\"95% CI\", color=\"C2\")\n",
    "ax.vlines(ci_upper, [0], 15, lw=2.5, linestyle=\"dotted\", color=\"C2\")\n",
    "\n",
    "ax.hist(\n",
    "    bootstrap_train_accuracies, bins=7, color=\"#0080ff\", edgecolor=\"none\", alpha=0.3\n",
    ")\n",
    "plt.legend(loc=\"upper left\")\n",
    "\n",
    "plt.xlim([0.8, 1.1])\n",
    "plt.xlabel('Accuracy')\n",
    "plt.ylabel('Count')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.grid()\n",
    "plt.savefig(\"matplotlib-figures/bootstrap-hist-oob-t.pdf\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- (We will compare the interval to the normal approximation interval in a later section.)\n",
    "- Again, let's add our CI values to the Python dictionary for a comparison study later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "results[\"Method 2.1: Bootstrap, 1-sample CI\"] = {\n",
    "    \"Test accuracy\": bootstrap_train_mean,\n",
    "    \"Lower 95% CI\": ci_lower,\n",
    "    \"Upper 95% CI\": ci_upper,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2) Bootstrap Confidence Intervals Using the Percentile Method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Now, while the approach outlined in the previous section seems pretty straightforward if our samples follow a normal distribution, there is a more robust and general approach for utilizing the bootstrap samples, namely, the percentile method (section 2, *Bootstrapping and Uncertainties* of my \"[Model Evaluation, Model Selection, and Algorithm Selection in Machine Learning\"](https://arxiv.org/abs/1811.12808) article).\n",
    "- Here, we pick our lower and upper confidence bounds as follows:\n",
    "\n",
    "  - $\\text{ACC}_{lower} = \\alpha_{1}th$ percentile of the $\\text{ACC}_\\text{boot}$ distribution\n",
    "  - $\\text{ACC}_{upper} = \\alpha_{2}th$ percentile of the $\\text{ACC}_{boot}$ distribution,\n",
    "\n",
    "- Here, $\\alpha_1 = \\alpha$ and $\\alpha_2 = 1 - \\alpha$, and $\\alpha$ is our degree of confidence to compute the $100 \\times (1 - 2 \\times \\alpha)$ confidence interval. For instance, to compute a 95% confidence interval, we pick $\\alpha = 0.025$ to obtain the 2.5th and 97.5th percentiles of the *b* bootstrap samples distribution as our upper and lower confidence bounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8695652173913043 1.0\n"
     ]
    }
   ],
   "source": [
    "ci_lower = np.percentile(bootstrap_train_accuracies, 2.5)\n",
    "ci_upper = np.percentile(bootstrap_train_accuracies, 97.5)\n",
    "\n",
    "print(ci_lower, ci_upper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Let's have a look how it looks like in the context of our "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEYCAYAAABRMYxdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAiWUlEQVR4nO3de5hddXno8e87uTC5ySQwxFxsE5p0EEMMGgIIlkG0UFFJVVRO0RSjaevBO8jFcsTLeUDB86DgoQZaHDgWS4EItEihkQ0UIjFgCMEACRJiYJqQhECGXEjI7/wxi5hkJsnOzF6zZ9b+fp5nnr3X77cu735dK7yu2y9SSkiSJBVJXbUDkCRJqjQLHEmSVDgWOJIkqXAscCRJUuFY4EiSpMLpX+0AytHQ0JAmTJhQ7TAK5dVXX2XIkCHVDqNQzGnlmdPKM6eVZ04r65FHHlmTUmrs7nr6RIEzcuRIFixYUO0wCqVUKtHc3FztMArFnFaeOa08c1p55rSyIuK5SqzHS1SSJKlwLHAkSVLhWOBIkqTC6RP34HRm69atrFy5ks2bN1c7lF6lvr6esWPHMmDAgGqHIklS1fTZAmflypUMGzaMcePGERHVDqdXSCmxdu1aVq5cyfjx46sdjiRJVdNnL1Ft3ryZgw46yOJmJxHBQQcd5FktSVLNy7XAiYgvR8QTEbE4Im6MiPqIGBER90TE0uxzeDfWX8lwC8GcSJKUY4ETEWOALwBTU0qTgH7AJ4DzgbkppYnA3GxakiSpYvK+RNUfGBQR/YHBwAvAaUBL1t8CTM85htxEBJ/85Cd3TG/bto3GxkY+8IEPVDEqSZKUW4GTUnoeuBxYAbQCL6eU7gZGppRas3lagUPyiiFvQ4YMYfHixWzatAmAe+65hzFjxlQ5KkmSlNtTVNm9NacB44H1wL9GxJn7sfwsYBZAY2MjpVJpl/4DDzyQDRs2VCrcLjvppJO4+eabmT59Otdffz0f/vCHeeihh9iwYQOvvvoq5557Lk888QSvv/46F1xwAaeeeirPPfccs2bNYuPGjQBcfvnlHH300TzwwANccsklHHTQQfz2t79lypQpXHvttft9X83mzZs75Gt3bW1t+5xH+8ec7r8HVm5lzabEwYOCd4/t+GoDc1p55rTyzGnvlOdj4u8Fnk0pvQgQEbcC7wJWRcSolFJrRIwCVne2cEppNjAboKmpKe0+zseSJUsYNmzYjul/XfB7bn5kZcWC/+g7x3L61Lfsc75PfepTfOtb3+L0009nyZIl/M3f/A3z589n2LBhXHLJJZx88snccMMNrF+/nmnTpvHBD36QQw89lF/+8pfU19ezdOlSzjjjDBYsWMDgwYNZtGgRTzzxBKNHj+a4445j0aJFHH/88fsVe319PUceeeRe53HslMozp/vv6h/P4+Fn13H0+BFc1Hxsh35zWnnmtPLMae+UZ4GzAjgmIgYDm4CTgAXAq8AM4NLs87ZKbGzlS5t4+Nl1lVgVAMccelBZ802ePJnly5dz44038v73v3+Xvrvvvpvbb7+dyy+/HGg/s7JixQpGjx7N2WefzcKFC+nXrx9PP/30jmWmTZvG2LFjAZgyZQrLly/f7wJHkqRal1uBk1J6OCJuBh4FtgG/of2MzFDgpoiYSXsRdHoltjd2+CCOHj+iEqvasb5yfehDH+Kcc86hVCqxdu3aHe0pJW655Raampp2mf/iiy9m5MiRPPbYY2zfvp36+vodfQcccMCO7/369WPbtm3d+BWSJNWmXN9knFL6BvCN3Zq30H42p6JOn/qWsi4p5eHTn/40Bx54IEccccQu12FPPvlkrrzySq688koigt/85jcceeSRvPzyy4wdO5a6ujpaWlp4/fXXqxK3JElF1WffZNybjB07li9+8Ysd2i+66CK2bt3K5MmTmTRpEhdddBEAn/vc52hpaeGYY47h6aefZsiQIT0dsiRJhdZnx6LqDdra2jq0NTc377jZbNCgQfz4xz/uMM/EiRNZtGjRjulLLrmkw7IAV111VWUDliSpRngGR5IkFY4FjiRJKhwLHEmSVDgWOJIkqXAscCRJUuFY4EiSpMKxwOmGH/zgB0yaNIm3ve1tXHHFFTvaL774YsaMGcOUKVOYMmUKd955JwAPPvggkydP5qijjmLZsmUArF+/npNPPpmUUqfb2Lp1K+effz4TJ05k0qRJTJs2jV/84hcAjBs3jjVr1uT7IyVJ6oN8D04XLV68mGuuuYb58+czcOBATjnlFE499VQmTpwIwJe//GXOOeecXZb5/ve/zy233MLy5cu5+uqr+f73v8+3v/1tLrzwwj2OGH7RRRfR2trK4sWLOeCAA1i1ahX33Xdf7r9PkqS+zDM4XbRkyRKOOeYYBg8eTP/+/TnhhBOYM2fOXpcZMGAAmzZtYuPGjQwYMIBnnnmG559/nhNOOKHT+Tdu3Mg111zDlVdeuWOMqpEjR/Kxj32s4r9HkqQiKUyB8/NlP+esu87irLvO6tB37n3nctZdZ3Ht49fu0v7kuid3LPPkuid36Vuzae+XfiZNmsT999/P2rVr2bhxI3feeSe///3vd/RfddVVTJ48mU9/+tO89NJLAFxwwQXMmjWLK664grPPPpuvf/3rfPvb397jNpYtW8Yf/dEf8aY3vWmfv1+SJP1BYQqcF9peYMGqBSxYtaBD36IXF7Fg1QKeffnZXdo3vLZhxzIbXtuwS9+W17fsdXtvfetbOe+883jf+97HKaecwtvf/nb692+/4vd3f/d3PPPMMyxcuJBRo0bx1a9+FYApU6bwq1/9invvvZff/e53jB49mpQSH//4xznzzDNZtWpVd1IgSZIyhSlwRg8dzdSRU5k6cmqHvsmNk5k6cirjDxy/S/uwgcN2LDNs4LBd+g7od8A+tzlz5kweffRR7r//fkaMGLHj/puRI0fSr18/6urq+OxnP8v8+fN3WS6lxHe+8x0uuugivvnNb/LNb36TM888kx/+8Ie7zDdhwgRWrFjBhg27Fl+SJGnvCnOT8fQJ05k+YXqnfZedcFmn7YeNOIzrTrmu076DBx28z22uXr2aQw45hBUrVnDrrbcyb948AFpbWxk1ahQAc+bMYdKkSbss19LSwqmnnsrw4cPZuHEjdXV11NXVsXHjxl3mGzx4MDNnzuQLX/gCP/7xjxk4cCCtra3MnTuXM888c5/xSZJUqwpT4FTDRz7yEdauXcuAAQP40Y9+xPDhwwH42te+xsKFC4kIxo0bt8uI4hs3bqSlpYW7774bgK985St85CMfYeDAgdx4440dtvGd73yHv//7v+fwww+nvr6eIUOG8K1vfatnfqAkSX2UBU43PPDAA52233DDDXtcZvDgwdx77707pt/97nfz+OOP73H+gQMH8r3vfY/vfe97HfqWL19efrCSJNWQwtyDI0mS9AYLHEmSVDi5FTgR0RQRC3f6eyUivhQRIyLinohYmn0O7+o29jS8QS0zJ5Ik5VjgpJSeSilNSSlNAd4JbATmAOcDc1NKE4G52fR+q6+vZ+3atf4HfScpJdauXUt9fX21Q5Ekqap66ibjk4BnUkrPRcRpQHPW3gKUgPP2d4Vjx45l5cqVvPjiixULsgjq6+sZO3ZstcOQJKmqeqrA+QTwxjPQI1NKrQAppdaIOKQrKxwwYADjx4/f94ySJKnmRN6XeCJiIPAC8LaU0qqIWJ9Satip/6WUUof7cCJiFjALoLGx8Z033XRTrnHWmra2NoYOHVrtMArFnO6/Sx7exFMvbadpeB0XHD2oQ785rTxzWnnmtLJOPPHER1JKHYcl2E89cQbnL4BHU0pvDLS0KiJGZWdvRgGrO1sopTQbmA3Q1NSUmpubeyDU2lEqlTCnlWVO99/VT82Dl9bR0NBAc/OxHfrNaeWZ08ozp71TTzwmfgZ/uDwFcDswI/s+A7itB2KQJEk1JNcCJyIGA+8Dbt2p+VLgfRGxNOu7NM8YJElS7cn1ElVKaSNw0G5ta2l/qkqSJCkXvslYkiQVjgWOJEkqHAscSZJUOBY4kiSpcCxwJElS4VjgSJKkwrHAkSRJhWOBI0mSCscCR5IkFY4FjiRJKhwLHEmSVDgWOJIkqXAscCRJUuFY4EiSpMKxwJEkSYVjgSNJkgrHAkeSJBWOBY4kSSocCxxJklQ4uRY4EdEQETdHxJMRsSQijo2IERFxT0QszT6H5xmDJEmqPf1zXv8PgLtSSh+NiIHAYOBCYG5K6dKIOB84Hzgv5zgkqUfc8VS1I9i7YdUOQOohuZ3BiYg3AX8G/CNASum1lNJ64DSgJZutBZieVwySJKk25XmJ6lDgReC6iPhNRFwbEUOAkSmlVoDs85AcY5AkSTUoUkr5rDhiKvAr4LiU0sMR8QPgFeDzKaWGneZ7KaXU4T6ciJgFzAJobGx850033ZRLnLWqra2NoUOHVjuMQjGn+++Shzfx1EvbaRpexwVHD+rQ3xdz+vKWakewd/229r2c9nZ9cT/tzU488cRHUkpTu7uePO/BWQmsTCk9nE3fTPv9NqsiYlRKqTUiRgGrO1s4pTQbmA3Q1NSUmpubcwy19pRKJcxpZZnT/Xf1U/PgpXU0NDTQ3Hxsh/6+mNPefg/O0Na+l9Peri/up7Ugt0tUKaX/Bn4fEU1Z00nAb4HbgRlZ2wzgtrxikCRJtSnvp6g+D/w0e4Lqd8BZtBdVN0XETGAFcHrOMUiSpBqTa4GTUloIdHYd7aQ8tytJkmqbbzKWJEmFY4EjSZIKxwJHkiQVjgWOJEkqHAscSZJUOBY4kiSpcCxwJElS4VjgSJKkwrHAkSRJhWOBI0mSCscCR5IkFY4FjiRJKhwLHEmSVDgWOJIkqXAscCRJUuFY4EiSpMKxwJEkSYVjgSNJkgrHAkeSJBVO/zxXHhHLgQ3A68C2lNLUiBgB/AswDlgOfCyl9FKecUiSpNrSE2dwTkwpTUkpTc2mzwfmppQmAnOzaUmSpIqpxiWq04CW7HsLML0KMUiSpALLu8BJwN0R8UhEzMraRqaUWgGyz0NyjkGSJNWYXO/BAY5LKb0QEYcA90TEk+UumBVEswAaGxsplUo5hVib2trazGmFmdP9t379puxzfae564s53b6l2hHsXdvWvpfT3q4v7qe1INcCJ6X0Qva5OiLmANOAVRExKqXUGhGjgNV7WHY2MBugqakpNTc35xlqzSmVSpjTyjKn++/qp+bBS+toaGigufnYDv19Mad3PFXtCPZuaGvfy2lv1xf301qQ2yWqiBgSEcPe+A78ObAYuB2Ykc02A7gtrxgkSVJtyvMMzkhgTkS8sZ1/TindFRG/Bm6KiJnACuD0HGOQJEk1KLcCJ6X0O+DtnbSvBU7Ka7uSJEm+yViSJBWOBY4kSSocCxxJklQ4FjiSJKlwLHAkSVLhlFXgRMRx5bRJkiT1BuWewbmyzDZJkqSq2+t7cCLiWOBdQGNEfGWnrjcB/fIMTJIkqav29aK/gcDQbL5hO7W/Anw0r6AkSZK6Y68FTkrpPuC+iPhJSum5HopJUpX11ICRazf94bOzbW7f0rH9g035xyWp7yt3qIYDImI2MG7nZVJK78kjKEmSpO4ot8D5V+AfgGuB1/MLR5IkqfvKLXC2pZSuzjUSSZKkCin3MfE7IuJzETEqIka88ZdrZJIkSV1U7hmcGdnnuTu1JeDQyoYjSZLUfWUVOCml8XkHIkmSVCllFTgR8anO2lNK11c2HEmSpO4r9xLVUTt9rwdOAh4FLHAkSVKvU+4lqs/vPB0RBwI35BKRJElSN5X7FNXuNgITKxmIJElSpZR7D84dtD81Be2DbL4VuKnMZfsBC4DnU0ofyB4v/xfa34q8HPhYSuml/QtbUq3qqWEkJPVt5d6Dc/lO37cBz6WUVpa57BeBJbSPQA5wPjA3pXRpRJyfTZ9X5rokSZL2qaxLVNmgm0/SPqL4cOC1cpaLiLHAqbQP8fCG04CW7HsLML3MWCVJkspS7iWqjwGXASUggCsj4tyU0s37WPQK4Gu0F0ZvGJlSagVIKbVGxCF72OYsYBZAY2MjpVKpnFBVpra2NnNaYUXK6fYtPbShzdlw4pvXs315qWP/lrbO29VlbVuLs5/2FkU69ouk3EtUXweOSimtBoiIRuA/gT0WOBHxAWB1SumRiGje38BSSrOB2QBNTU2puXm/V6G9KJVKmNPKKlJOe+w+lyXzgHVQ30DduGM7dG9fXqJuXHMPBVMbhrYWZz/tLYp07BdJuQVO3RvFTWYt+768dRzwoYh4P+3vznlTRPw/YFVEjMrO3owCVu91LZIkSfup3MfE74qI/4iIv46Ivwb+HbhzbwuklC5IKY1NKY0DPgH8MqV0JnA7fxjbagZwW5cilyRJ2oO9nsGJiAm03zNzbkR8GDie9ntw5gE/7eI2LwVuioiZwArg9C6uR5IkqVP7ukR1BXAhQErpVuBWgIiYmvV9sJyNpJRKtN+gTEppLe1DPUiSJOViX5eoxqWUFu3emFJaQPuL+iRJknqdfRU49XvpG1TJQCRJkiplXwXOryPis7s3ZvfPPJJPSJIkSd2zr3twvgTMiYi/4g8FzVRgIPCXOcYlSZLUZXstcFJKq4B3RcSJwKSs+d9TSr/MPTJJkqQuKutFfymle4F7c45FkiSpIsp90Z8kSVKfYYEjSZIKxwJHkiQVjgWOJEkqHAscSZJUOBY4kiSpcCxwJElS4VjgSJKkwrHAkSRJhWOBI0mSCscCR5IkFY4FjiRJKhwLHEmSVDi5FTgRUR8R8yPisYh4IiK+mbWPiIh7ImJp9jk8rxgkSVJtyvMMzhbgPSmltwNTgFMi4hjgfGBuSmkiMDebliRJqpjcCpzUri2bHJD9JeA0oCVrbwGm5xWDJEmqTf3zXHlE9AMeASYAP0opPRwRI1NKrQAppdaIOGQPy84CZgE0NjZSKpXyDLXmtLW1mdMKK1JOt2/poQ1t3pR9rmf78lLH/i1tnbery9q2Fmc/7S2KdOwXSa4FTkrpdWBKRDQAcyJi0n4sOxuYDdDU1JSam5tzibFWlUolzGllFSmndzzVQxtaMg9YB/UN1I07tkP39uUl6sY191AwtWFoa3H2096iSMd+kfTIU1QppfVACTgFWBURowCyz9U9EYMkSaodeT5F1ZiduSEiBgHvBZ4EbgdmZLPNAG7LKwZJklSb8rxENQpoye7DqQNuSin9W0TMA26KiJnACuD0HGOQJEk1KLcCJ6W0CDiyk/a1wEl5bVfqC3rsHhdJqlG+yViSJBWOBY4kSSocCxxJklQ4FjiSJKlwLHAkSVLhWOBIkqTCscCRJEmFY4EjSZIKxwJHkiQVjgWOJEkqHAscSZJUOBY4kiSpcCxwJElS4VjgSJKkwrHAkSRJhWOBI0mSCscCR5IkFY4FjiRJKhwLHEmSVDi5FTgR8ZaIuDcilkTEExHxxax9RETcExFLs8/hecUgSZJqU55ncLYBX00pvRU4BvifEXE4cD4wN6U0EZibTUuSJFVMbgVOSqk1pfRo9n0DsAQYA5wGtGSztQDT84pBkiTVpv49sZGIGAccCTwMjEwptUJ7ERQRh+xhmVnALIDGxkZKpVJPhFoz2trazGmF7U9Ot2/JN5Y+Y/Om7HM925eXOvZvaeu8XV3WttVjv9L897R3yr3AiYihwC3Al1JKr0REWcullGYDswGamppSc3NzbjHWolKphDmtrP3J6R1P5RtLn7FkHrAO6huoG3dsh+7ty0vUjWvu8bCKbGirx36l+e9p75TrU1QRMYD24uanKaVbs+ZVETEq6x8FrM4zBkmSVHvyfIoqgH8ElqSU/s9OXbcDM7LvM4Db8opBkiTVpjwvUR0HfBJ4PCIWZm0XApcCN0XETGAFcHqOMUiSpBqUW4GTUvovYE833JyU13YlSZJ8k7EkSSocCxxJklQ4PfIeHKkaevpR7O1bfPxbknoLz+BIkqTCscCRJEmFY4EjSZIKxwJHkiQVjgWOJEkqHJ+iUlX8fNnPeaHtBUYPHc30CdOrHY5UEx5Z/XM2rn+I9cvWe9yp8DyDo6q4bdltXP3Y1dy2zKHIpJ7y6Iu38YuXf+Fxp5rgGRxVxWEjDtvlU1L+Rg05jH6vrfe4U02wwFFVnDftvGqHINWcD4w7j2EHlGie1lztUKTceYlKkiQVjgWOJEkqHAscVcWaTWt4vu151mxaU+1QpJqx4bU1rN221uNONcF7cFQV5953LgtWLWDqyKlcd8p11Q5Hqgk/W3ouz25YwNT7PO5UfJ7BkSRJheMZHFXFzCNmctqE0zh40MHVDkWqGSeMmcm7Xmzi+COOr3YoUu4scFQVx4/xH1ipp/1pw/EM27TN4081IbdLVBHxTxGxOiIW79Q2IiLuiYil2efwvLYvSZJqV5734PwEOGW3tvOBuSmlicDcbFqSJKmicrtElVK6PyLG7dZ8GtCcfW8BSoCvtK1BDrYp9TwH21Qt6el7cEamlFoBUkqtEXHInmaMiFnALIDGxkZKpVLPRFgj2traqprTlv9uYdmWZUw4YAINKxty2cb2Lbmsds+2tLF9eamHN9rHbd6Ufa7vPHfmtKIeWdvC8q3LWLpgaW7HXS2q9r+n6lyvvck4pTQbmA3Q1NSUmpubqxtQwZRKJaqZ05a7WmAVNDQ05BbHHU/lsto92r68RN245p7daF+3ZB6wDuobqBt3bIduc1pZ8WoLbM33uKtF1f73VJ3r6QJnVUSMys7ejAJW9/D21Uv4kjGp5332bdcxrNX/GKs29PSL/m4HZmTfZwC39fD2JUlSDcjzMfEbgXlAU0SsjIiZwKXA+yJiKfC+bFqSJKmi8nyK6ow9dJ2U1zYlSZKgF99krGI7975zWfTiIiY3TuayEy6rdjhSTbhx6bk8v34+0+6b5nGnwnOwTVXFmk1reOHVF1izaU21Q5FqRttra1j3+jqPO9UEz+CoKo4bcxyjh45m/IHjqx2KVDP+tOE4GrfXceyYjo/kS0VjgaOq+MwRn6l2CFLNOWHMZxhWN4HmI5qrHYqUOy9RSZKkwrHAkSRJheMlKlXFk+ueZMNrGxg2cBiHjTis2uFINeGFV5+kbvNS3rzuzR53KjwLHFXFd+d/lwWrFjB15FSHbZB6yL8v/y7PbljAQ/Mf8rhT4VngqMu6M5jl2o1/+OzpQTElScVngaOqOHXceWzetoH6/sOqHYpUM04ddx51qx/g3dPeXe1QpNxZ4KgqRg/x+r/U00YPOYxh9f/t/TeqCT5FJUmSCscCR5IkFY6XqHqxPG++3b6lujf33vf8tby46VkaB43nhDG+1VjqCfc9fy0vrZvHsseX+TZxFZ5ncFQVT69/kEfX3M7T6x+sdihSzXh6/YPMf3U+Dz7vcafis8BRVQwdeDANA0czdODB1Q5FqhlDBx7MiH4jOHiQx52Kz0tUqoozJl5W7RCkmnPGxMsY1lqi+YTmaoci5c4zOJIkqXAscCRJUuFUpcCJiFMi4qmIWBYR51cjBkmSVFw9XuBERD/gR8BfAIcDZ0TE4T0dh6rrmifO4sJfHcE1T5xV7VCkmnHNE2fx+ec+z1l3edyp+KpxBmcasCyl9LuU0mvAz4DTqhCHJEkqqEgp9ewGIz4KnJJS+kw2/Ung6JTS2bvNNwuYlU1OAhb3aKDFdzCwptpBFIw5rTxzWnnmtPLMaWU1pZS6PRJzNR4Tj07aOlRZKaXZwGyAiFiQUpqad2C1xJxWnjmtPHNaeea08sxpZUXEgkqspxqXqFYCb9lpeizwQhXikCRJBVWNAufXwMSIGB8RA4FPALdXIQ5JklRQPX6JKqW0LSLOBv4D6Af8U0rpiX0sNjv/yGqOOa08c1p55rTyzGnlmdPKqkg+e/wmY0mSpLz5JmNJklQ4FjiSJKlwqlrg7GvIhog4MCLuiIjHIuKJiDir3GVrVTdzujwiHo+IhZV6TK8Iysjp8IiYExGLImJ+REwqd9la1c2cup92IiL+KSJWR0Sn7wyLdj/Mcr4oIt6xU5/76W66mU/30U6UkdPDImJeRGyJiHN269v/fTSlVJU/2m8wfgY4FBgIPAYcvts8FwLfzb43Auuyefe5bC3+dSen2fRy4OBq/47e9FdmTi8DvpF9PwyYW+6ytfjXnZxm0+6nnef1z4B3AIv30P9+4Be0v4vsGODhcv/3qMW/ruYz63Mf7VpODwGOAv43cM5O7V3aR6t5BqecIRsSMCwiAhhK+3+Mt5W5bC3qTk7VuXJyejgwFyCl9CQwLiJGlrlsLepOTrUHKaX7aT+e9+Q04PrU7ldAQ0SMwv20U93Ip/ZgXzlNKa1OKf0a2LpbV5f20WoWOGOA3+80vTJr29lVwFtpfxHg48AXU0rby1y2FnUnp9Be/NwdEY9kQ2WovJw+BnwYICKmAX9M+wss3U87152cgvtpV+0p7+6nXbO3vLmPVlaX9tFqDNXwhnKGbDgZWAi8B/gT4J6IeKDMZWtRl3OaUnoFOC6l9EJEHJK1P5lV3LWsnJxeCvwgIhbSXjT+hvazYu6nnetOTsH9tKv2lHf3067ZW97cRyurS/toNc/glDNkw1nArdkpwGXAs7Rfj3e4h851J6eklF7IPlcDc2g/LVjr9pnTlNIrKaWzUkpTgE/Rfm/Ts+UsW6O6k1P3067bU97dT7tmj3lzH624Lu2j1SxwyhmyYQVwEkB2/b0J+F2Zy9aiLuc0IoZExLCsfQjw5ziCO5SR04hoyPoAPgPcn50Rcz/tXJdz6n7aLbcDn8qe/jkGeDml1Ir7aVd1mk/30Vx0aR+t2iWqtIchGyLib7P+fwC+DfwkIh6n/RTVeSmlNQCdLVuN39GbdCenEXEoMKf93mP6A/+cUrqrKj+kFykzp28Fro+I14HfAjP3tmw1fkdv0p2cAiNxP+1URNwINAMHR8RK4BvAANiR0ztpf/JnGbCR9rO57qd70NV84j66R/vKaUS8GVgAvAnYHhFfov1pqVe6so86VIMkSSoc32QsSZIKxwJHkiQVjgWOJEkqHAscSZJUOBY4kiSpcCxwJFVERPxlRKSIOKzasUiSBY6kSjkD+C/aX8KVi4jol9e6JRWLBY6kbouIocBxtL+Q7xNZW7+IuDwiHo+IRRHx+az9qIh4KCIei4j5ETEsIv46Iq7aaX3/FhHN2fe2iPhWRDwMHBsR/ysifh0RiyNidmRvVIuICRHxn9l6H42IP4mIGyLitJ3W+9OI+FBP5UVS9VjgSKqE6cBdKaWngXUR8Q5gFjAeODKlNBn4afaa9X+hfRT7twPvBTbtY91DgMUppaNTSv8FXJVSOiqlNAkYBHwgm++nwI+y9b4LaAWuJXvDbEQcmLXfWakfLan3ssCRVAlnAD/Lvv8sm34v8A8ppW0AKaV1tI991ppS+nXW9sob/XvxOnDLTtMnRsTD2XAj7wHelo39MyalNCdb7+aU0saU0n3AhGxU5zOAW8rYnqQCqNpYVJKKISIOor3QmBQRifaxYhLwSPa5y+ydtAFsY9f/w1W/0/fNKaXXs23VA/8XmJpS+n1EXJzNG3sJ8Qbgr2i/dPbpMn+WpD7OMziSuuujwPUppT9OKY1LKb0FeBZ4FPjbiOgPEBEjgCeB0RFxVNY2LOtfDkyJiLqIeAswbQ/beqPwWZPd9/NRaD8TBKyMiOnZeg+IiMHZvD8BvpTNV/ODSEq1wgJHUnedAczZre0WYDSwAlgUEY8B/yOl9BrwceDKrO0e2ouWB2kvih4HLqe9OOogpbQeuCab7+fAr3fq/iTwhYhYBDwEvDlbZhWwBLium79TUh/iaOKSCi07k/M48I6U0svVjkdSz/AMjqTCioj30n5Z7EqLG6m2eAZHkiQVjmdwJElS4VjgSJKkwrHAkSRJhWOBI0mSCscCR5IkFc7/BwGRJ+nAVITvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "ax.vlines(bootstrap_train_mean, [0], 80, lw=2.5, linestyle=\"-\", label=\"Mean\")\n",
    "\n",
    "ax.vlines(ci_lower, [0], 15, lw=2.5, linestyle=\"dotted\", label=\"95% CI\", color=\"C2\")\n",
    "ax.vlines(ci_upper, [0], 15, lw=2.5, linestyle=\"dotted\", color=\"C2\")\n",
    "\n",
    "ax.hist(\n",
    "    bootstrap_train_accuracies, bins=7, color=\"#0080ff\", edgecolor=\"none\", alpha=0.3\n",
    ")\n",
    "plt.legend(loc=\"upper left\")\n",
    "\n",
    "plt.xlim([0.8, 1.1])\n",
    "plt.xlabel('Accuracy')\n",
    "plt.ylabel('Count')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.grid()\n",
    "plt.savefig(\"matplotlib-figures/bootstrap-hist-oob-percentile.pdf\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "results[\"Method 2.2: Bootstrap, percentile\"] = {\n",
    "    \"Test accuracy\": bootstrap_train_mean,\n",
    "    \"Lower 95% CI\": ci_lower,\n",
    "    \"Upper 95% CI\": ci_upper,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3) Reweighting the Boostrap Samples via the .632 Bootstrap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Building on the previously introduced bootstrap percentile method, a slightly improved version is the [.632 Bootstrap](https://scholar.google.com/scholar?hl=en&as_sdt=0%2C50&q=Estimating+the+error+rate+of+a+prediction+rule%3A+improvement+on+cross-validation&btnG=).\n",
    "\n",
    "- Skipping over the technical details the previously introduced bootstrap method has a small pessimistic bias, which means that it reports a test accuracy that is slightly worse than the true generalization accuracy of the model; the so-called .632 bootstrap aims to address that.\n",
    "\n",
    "- For a more detailed discussion please see section 2, *Bootstrapping and Uncertainties* of my \"[Model Evaluation, Model Selection, and Algorithm Selection in Machine Learning\"](https://arxiv.org/abs/1811.12808) article).\n",
    "\n",
    "- We are skipping the formulas and jump directly into the code implementation (in practice, I recommend using my implementation in [mlxtend](http://rasbt.github.io/mlxtend/user_guide/evaluate/bootstrap_point632_score/)).\n",
    "\n",
    "- In a nutshell, you can think of it as a reweighted version of the bootstrap method we used earlier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9677367193053941"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rng = np.random.RandomState(seed=12345)\n",
    "idx = np.arange(y_train.shape[0])\n",
    "\n",
    "bootstrap_train_accuracies = []\n",
    "bootstrap_rounds = 200\n",
    "weight = 0.632\n",
    "\n",
    "for i in range(bootstrap_rounds):\n",
    "\n",
    "    train_idx = rng.choice(idx, size=idx.shape[0], replace=True)\n",
    "    valid_idx = np.setdiff1d(idx, train_idx, assume_unique=False)\n",
    "\n",
    "    boot_train_X, boot_train_y = X_train[train_idx], y_train[train_idx]\n",
    "    boot_valid_X, boot_valid_y = X_train[valid_idx], y_train[valid_idx]\n",
    "\n",
    "    clf.fit(boot_train_X, boot_train_y)\n",
    "    valid_acc = clf.score(boot_valid_X, boot_valid_y)\n",
    "    # predict training accuracy on the whole training set\n",
    "    # as ib the original .632 boostrap paper\n",
    "    # in Eq (6.12) in\n",
    "    #    \"Estimating the Error Rate of a Prediction Rule: Improvement\n",
    "    #     on Cross-Validation\"\n",
    "    #     by B. Efron, 1983, https://doi.org/10.2307/2288636\n",
    "    train_acc = clf.score(X_train, y_train)\n",
    "\n",
    "    acc = weight * train_acc + (1.0 - weight) * valid_acc\n",
    "\n",
    "    bootstrap_train_accuracies.append(acc)\n",
    "\n",
    "bootstrap_train_mean = np.mean(bootstrap_train_accuracies)\n",
    "bootstrap_train_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9221417322834646 1.0\n"
     ]
    }
   ],
   "source": [
    "ci_lower = np.percentile(bootstrap_train_accuracies, 2.5)\n",
    "ci_upper = np.percentile(bootstrap_train_accuracies, 97.5)\n",
    "\n",
    "print(ci_lower, ci_upper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "results[\"Method 2.3: Bootstrap, .632\"] = {\n",
    "    \"Test accuracy\": bootstrap_train_mean,\n",
    "    \"Lower 95% CI\": ci_lower,\n",
    "    \"Upper 95% CI\": ci_upper,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4) Taking the Reweighting One Step Further: The .632+ Bootstrap "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [.632+ Bootstrap](https://scholar.google.com/scholar_lookup?&title=Improvements%20on%20Cross-Validation%3A%20The%20.632%2B%20Bootstrap%20Method&journal=J%20Am%20Stat%20Assoc&doi=10.1080%2F01621459.1997.10474007&volume=92&issue=438&pages=548-560&publication_year=1997&author=Efron%2CB&author=Tibshirani%2CR) is an improvement over the .632 Bootstrap we implemented above. In a nutshell, the main difference is that the weighting terms are computed from the so-called no-information rate rather than being fixed.\n",
    "\n",
    "- Again, we are skipping the formulas and jump directly into the code implementation (in practice, I recommend using my implementation in [mlxtend](http://rasbt.github.io/mlxtend/user_guide/evaluate/bootstrap_point632_score/)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9683445668115584"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from itertools import product\n",
    "\n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "\n",
    "def no_information_rate(targets, predictions, loss_fn):\n",
    "    combinations = np.array(list(product(targets, predictions)))\n",
    "    return loss_fn(combinations[:, 0], combinations[:, 1])\n",
    "\n",
    "\n",
    "rng = np.random.RandomState(seed=12345)\n",
    "idx = np.arange(y_train.shape[0])\n",
    "\n",
    "bootstrap_train_accuracies = []\n",
    "bootstrap_rounds = 200\n",
    "weight = 0.632\n",
    "\n",
    "\n",
    "for i in range(bootstrap_rounds):\n",
    "\n",
    "    train_idx = rng.choice(idx, size=idx.shape[0], replace=True)\n",
    "    valid_idx = np.setdiff1d(idx, train_idx, assume_unique=False)\n",
    "\n",
    "    boot_train_X, boot_train_y = X_train[train_idx], y_train[train_idx]\n",
    "    boot_valid_X, boot_valid_y = X_train[valid_idx], y_train[valid_idx]\n",
    "\n",
    "    clf.fit(boot_train_X, boot_train_y)\n",
    "    train_acc = clf.score(X_train, y_train)\n",
    "    valid_acc = clf.score(boot_valid_X, boot_valid_y)\n",
    "\n",
    "    gamma = no_information_rate(y, clf.predict(X), accuracy_score)\n",
    "    R = (valid_acc - train_acc) / (gamma - train_acc)\n",
    "\n",
    "    weight = 0.632 / (1 - 0.368 * R)\n",
    "\n",
    "    acc = weight * train_acc + (1.0 - weight) * valid_acc\n",
    "\n",
    "    bootstrap_train_accuracies.append(acc)\n",
    "\n",
    "bootstrap_train_mean = np.mean(bootstrap_train_accuracies)\n",
    "bootstrap_train_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9248753660725227 1.0\n"
     ]
    }
   ],
   "source": [
    "ci_lower = np.percentile(bootstrap_train_accuracies, 2.5)\n",
    "ci_upper = np.percentile(bootstrap_train_accuracies, 97.5)\n",
    "\n",
    "print(ci_lower, ci_upper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "results[\"Method 2.4: Bootstrap, .632+\"] = {\n",
    "    \"Test accuracy\": bootstrap_train_mean,\n",
    "    \"Lower 95% CI\": ci_lower,\n",
    "    \"Upper 95% CI\": ci_upper,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3) Bootstrapping the Test Set Predictions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Next, let's look at another way we can construct confidence intervals involving bootstrapping, namely, bootstrapping the test set.\n",
    "- I saw this method being used in the article [Machine Learning for Scent: Learning Generalizable Perceptual Representations of Small Molecules](https://arxiv.org/abs/1910.10685).\n",
    "- Here, in contrast to the other bootstrap methods we covered previously, we keep the model fixed and just resample the test set (instead of the training set). In a nutshell, this method avoid retraining the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9597826086956522"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.fit(X_train, y_train)\n",
    "\n",
    "predictions_test = clf.predict(X_test)\n",
    "acc_test = np.mean(predictions_test == y_test)\n",
    "\n",
    "rng = np.random.RandomState(seed=12345)\n",
    "idx = np.arange(y_test.shape[0])\n",
    "\n",
    "test_accuracies = []\n",
    "\n",
    "for i in range(200):\n",
    "\n",
    "    pred_idx = rng.choice(idx, size=idx.shape[0], replace=True)\n",
    "    acc_test_boot = np.mean(predictions_test[pred_idx] == y_test[pred_idx])\n",
    "    test_accuracies.append(acc_test_boot)\n",
    "\n",
    "bootstrap_train_mean = np.mean(test_accuracies)\n",
    "bootstrap_train_mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- We can then use the familiar percentile approach to get the confidence interval:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8260869565217391 1.0\n"
     ]
    }
   ],
   "source": [
    "ci_lower = np.percentile(test_accuracies, 2.5)\n",
    "ci_upper = np.percentile(test_accuracies, 97.5)\n",
    "\n",
    "print(ci_lower, ci_upper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "results[\"Method 3: Bootstrap test set\"] = {\n",
    "    \"Test accuracy\": bootstrap_train_mean,\n",
    "    \"Lower 95% CI\": ci_lower,\n",
    "    \"Upper 95% CI\": ci_upper,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4) Confidence Intervals from Retraining Models with Different Random Seeds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- In deep learning, it is quite common to retrain a model with different random seeds. How can we construct confidence interval from these experiments?\n",
    "\n",
    "- Assuming that the sample means are normal distributed, we can adopt the approach from earlier:\n",
    "\n",
    "$$\n",
    "\\bar{x} \\pm z \\times \\text{SE}.\n",
    "$$\n",
    "\n",
    "- And since we usually deal with a smaller dataset size, a t-distribution is more appropriate, similar to what we used in the first bootstrap section (before we introduced the percentile method)\n",
    "\n",
    "- Also, if we are interested in the average accruacy, $\\overline{ACC}_{\\text{test}}$, we can technically make the argument that each $\\text{ACC}_{\\text{test}}$ corresponding to a different random seed is a sample, and the number of random seeds we evaluate would be the sample size $n$ so that we compute\n",
    "\n",
    "$$\n",
    "\\overline{ACC}_{\\text{test}} \\pm t \\times \\text{SE},\n",
    "$$\n",
    "\n",
    "with\n",
    "\n",
    "$$\n",
    "\\text{SE} = \\frac{\\text{SD}}{\\sqrt{n}}.\n",
    "$$\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Note that for a decision tree classifier, the random seed usually doesn't really make difference in practice, so the following experiment is not very interesting. However, I am including the code just for the sake of completeness so that you can get an idea of how it works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9565217391304348"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_accuracies = []\n",
    "rounds = 5\n",
    "\n",
    "\n",
    "for i in range(rounds):\n",
    "\n",
    "    clf = DecisionTreeClassifier(random_state=i)\n",
    "\n",
    "    clf.fit(X_train, y_train)\n",
    "    acc = clf.score(X_test, y_test)\n",
    "    test_accuracies.append(acc)\n",
    "\n",
    "test_mean = np.mean(test_accuracies)\n",
    "test_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9565217391304348 0.9565217391304348\n"
     ]
    }
   ],
   "source": [
    "confidence = 0.95  # Change to your desired confidence level\n",
    "t_value = scipy.stats.t.ppf((1 + confidence) / 2.0, df=rounds - 1)\n",
    "\n",
    "sd = np.std(test_accuracies, ddof=1)\n",
    "se = sd / np.sqrt(rounds)\n",
    "\n",
    "ci_length = t_value * se\n",
    "\n",
    "ci_lower = test_mean - ci_length\n",
    "ci_upper = test_mean + ci_length\n",
    "\n",
    "print(ci_lower, ci_upper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- As suspected, the test accuracies are all identical. However, in the context of training deep neural networks, this is another viable method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Comparing the Different Confidence Interval Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Given that there are so many methods out there, which one should we use? \n",
    "\n",
    "- It's a tricky to give a general recommendation since there are two aspects to it: practicality and accuracy.\n",
    "\n",
    "- Let's look at the practicality first.\n",
    "\n",
    "**Practicality**\n",
    "\n",
    "- The normal approximation approach is great if we want a computationally cheap method for confidence intervals that avoids retraining the model in contrast to the bootstrap methods.\n",
    "\n",
    "- Similar to the normal approximation approach, the \"bootstrapping the test set\" method also avoid retraining the model, however, it requires that we have access to the model's test set predictions. In contrast, the normal approximation intervals can be computed just from the tabulated test set scores (and sizes) listed in a paper without rerunning additional experiments.\n",
    "\n",
    "- The other bootstrap methods are much more expensive because they involve retraining models on the training folds. Since a minimum of 200 bootstrap rounds is commonly recommended, this can be very expensive for bigger dataset and deep neural networks. Also, we do not get a single model in the end that we evaluate; to circumvent that, we could train the model on the training set and then evaluate it's performance by fitting 200 models on the bootsrap samples and computing the average performance. This works well for most machine learning classifiers, but we have to take care in case of deep learning models as they may not always converge, and the non-converging models can then produce misleading accuracy estimates.\n",
    "\n",
    "- The .632+ bootstrap might be the most accurate bootstrap method but it is computationally very expensive for large datasets. In the current implementation, it is likely not feasible for more than a few hundred training examples. Hence, if bootstrapping is used, the next best approach, the .632 bootstrap, might be a better alternative.\n",
    "\n",
    "- Computing the confidence intervals from different random seeds is another great option, however, it is only really useful for deep learning models. It's more expensive than Normal approximation and bootstrapping the test set, since it involves retraining the model. On the other hand, the results from different random seeds give us a good idea of the stability of the model. In fact, this is great for model comparisons if you do this for two different types of models. In this case, you can apply the following formula assuming unequal variances:\n",
    "\n",
    "\n",
    "$$\n",
    "\\left(\\overline{ACC}_{\\text{m1}} -\\overline{ACC}_{\\text{m2}} \\right) \\pm t \\sqrt{\\frac{\\text{SD}_{\\text{m1}}^{2}}{n_{\\text{m1}}}+\\frac{\\text{SD}_{\\text{m2}}^{2}}{n_{\\text{m2}}}},\n",
    "$$\n",
    "\n",
    "where $\\text{m1}$ and $\\text{m2}$ refer to model 1 and model 2, respectively. If the confidence interval does not contain 0, then the performance of the models is statistically significant.\n",
    "\n",
    "- So, from a practicality standpoint, we can rank the methods as follows, from most to least practical\n",
    "\n",
    "1. Normal approximation;\n",
    "2. Bootstrapping the test set;\n",
    "3. Confidence intervals for different random seeds (deep learning only);\n",
    "4. Bootstrapping training sets with the percentile method;\n",
    "5. Bootstrapping training sets with the t-interval (tied with 4.);\n",
    "6. .632 bootstrap (tied with 4.);\n",
    "7. .632+ bootstrap."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Next, let's look at the methods side by side:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['Method 1: Normal approximation', 'Method 2.1: Bootstrap, 1-sample CI', 'Method 2.2: Bootstrap, percentile', 'Method 2.3: Bootstrap, .632', 'Method 2.4: Bootstrap, .632+', 'Method 3: Bootstrap test set'])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAADQCAYAAAD4dzNkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5KElEQVR4nO3de5xVVf3/8ddbQBkuihfyYZhimGMGNCR5CdTBTMoLglbeSvGSWWqZecHsa6QVmH7NS99Q80eYeSdEFA0vMGiooDAw42284mWy0HRUYEQcPr8/1jq6OZwz58yNmT3zeT4e8/Dstfdee33WHuez19qbs2VmOOeccy5dNmnvBjjnnHOu6TyBO+eccynkCdw555xLIU/gzjnnXAp5AnfOOedSyBO4c845l0KewJ1znZKkEZJekLRS0lhJ90k6Ps+2AyWZpO4bu535xHZ/vr3bUYikcklvtHc7uiJP4M45ACR9UdJcSe9JelHSuMS6TIJbmfj5n8T6YyS9KekVSeWJ8kGSHpXUbeNGA8BFwB/NrI+ZzTSzb5nZDe3QjmaJ7X65mG3judm5rdvkOpYOc7XpnGs/ceR5F3AN8A1gP+BuScPM7PnEpv3M7OMc+04GvgLsDvwRGBxXXwWcZWYNbRxCLjsCT7fDcVNFUvfsc+rSwUfgzjmAXYHPAn8wswYzmwssAL5fxL5bA7Vm9ibwIPB5AEnfjuWPF6pA0g8kPSvpA0nPSPpKLP+ipApJdZKeljQmsc80Sf8naXbcb6GkQXHdS7Edd8fZgs1iPSfH9d0kXSbpbUkvAwdntWcLSf8vzirUSvpNZhZB0nhJ/4z7vxtnHb6V2HcrSX+R9K+4fmZi3SGSlsZ4HpU0tJE++WRUXSDWh+Muy2KsRxY6lqTlks6TVAWskvRLSdOzjn+lpKvi5xMS5+dlST9spN3nxT77QFKNpK/n29a1kJn5j//4Txf/AYYAKwElyh4A7oyfBwIG1AJvAH8BtonrNgGeB7YHDgWeAPoAS4Gtizj2d2K9XwUE7EwYPfcAXgR+AWwK7A98AJTG/aYB7wB7EGYTbwJuTdS7HDggsVwBnBw/nwo8B3wO2AqYF+PrHtfPBK4FegOfARYBP4zrxgNrgR8A3YAfAf/K9B0wG7gN2DLGsF8s/wqwAtgz7nd8bONmefrFgJ2LjPWTbYs5Vvy8NMZfEvt7NbB5XN8NeBPYKy4fDAyK52e/uO1X4rpy4I34uRR4Hfhs4vdmUHv/fnfWHx+BO+cgJLMVwDmSekg6kPCHuldc/zYhwe5ImCbvS0gimNk6QhKbDpxNSGwXAVcDQyTNkzRH0mByOxn4vZk9YcGLZvYqsBfhQmCymX1kYVbgHuDoxL4zzGyRhSngm4CyIuP9LnCFmb1uZu8AkzIrJG0LfAs408xWmdkK4A/AUYn9XzWzP1u4NXADsB2wraTt4r6nmtm7ZrbWzObHfX4AXGtmCy3MctwArIlxFqMpsRZzrKti/PWxv5cAY+O6/YHVFmdPzGy2mb0Uz8984H5gnxzHbQA2A3aT1MPMlpvZS0XG55rIE7hzDjNbS/jjfTDwb+DnwO2E0TZmttLMnjSzj83sP8DpwIGSNo/rHzKzvcxsP2AdMJwwaryRMGK9GLg+z+E/B+T6I/9Z4PV4gZDxKjAgsfzvxOfVhIRfjM8SRorJejMyo/834/RzHWE0/plcxzWz1fFjH0Is75jZuzmOuSPw80ydsd7PxbYUoymxFnOs17P2uZlPL46OicsASPqWpMclvRPrOgjYJvugZvYicCYwEVgh6VZJxcbnmsgTuHMOADOrMrP9zGxrMxtNuIe8KN/m8b9KFkoS4SG2nxD+wHeLo7sngHz3e18nTM9m+xfwOUnJv1M7EKbbW+pNQkJL1ptszxrCLYJ+8WdzM/tSEfW+DmwlqV+edb9N1NnPzHqZ2S3NDaJAOwodK/tVlHcA5ZK2B8YRE7ikzYC/A5cB25pZP+Bess79J5Wa3WxmIwkXEQZc0opxuQRP4M45ACQNldRTUi9JZxOmhafFdXtKKpW0iaStCU+XV5jZe1nVnAxUmtlS4L9AiaTdgFFAvn8SdT1wtqTdFewsaUdgIbAKODdO65cT7rHf2grh3g78RNL2krYEJmRWWHgY737gfyVtHmMeJGm/QpXGfe8D/iRpy9jufePqPwOnxr6UpN6SDpbUtxXi+Q/x4cHmHsvM3iI8J/AX4BUzezau2pQwLf4W8HF8YO/AXHXE35H9Y9L/EKgnTKu7NuAJ3DmX8X3CyHQF8HXgG2a2Jq77PPAPwkNkTxFGqMl70UjaBvgp8D8A8V7t6cBcwj9POyPXQc3sDuC3hBHfB4QHyLYys4+AMYR7ym8DfwKOM7PnWiHWPwNzgGWEe78zstYfR0hczwDvEu7vb1dk3d8nPOSWea7gTAAze5Jwb/qPsc4XCbcXWsNE4IY4Xf7dFhzrZuAAEtPnZvYBYUbl9ljXMcCsPPtvRvgnhW8Tpvw/Q3gI0bWBzFOTzjnnnEsRH4E755xzKeQJ3DnnnEshT+DOOedcCnkCd84551LIX2biuqx+/frZzjt37hc4rVq1it69e7d3M9pcV4izK8QIXSPOxYsXv21m/Vtajydw12Vtu+22PPnkk+3djDZVUVFBeXl5ezejzXWFOLtCjNA14pT0auGtCvMpdOeccy6FPIE755xzKZSqBB7fj3tjYrm7pLck3VNgvzJJByWWJ8avimxuO3LuL+lUSdXxHbz/jF8hWaiuivjO3KXxfbuntKBdBb/xSNKZknoV2q6Zx1+vn5u4bz9JP26ldoz3Fyg413nNrKxlxOS57DRhNiMmz2VmZWt8PX76pCqBE74XebCkkrj8DYp7sUEZ4e05be1mMxtiZmXA74HLi9zv2LjPCOASSZs28/jFfGXhmXz6isj1SOrWzONmlNH8fu4HtEoCJ3xlpCdw5zqhmZW1nD+jmtq6+vCC+rp6zp9R3SWTeBofYruP8MrD6YTvYr6F+F5aSb2J7yAmxDYxbn8R4aUKI/n0vb+7SaogvIXoCjO7KtZxFnBi3OZ6M7sill9A+H7k1wlf6r84u2Fm9n5isTcbvu2nkD6Ei5SGeMyjCUlZwGwzOy9fuaTJMcalwNPAKYTvLt4e6EZ4neO2hMQ2T9LbZjZK0krChcZowusH9ye8MKIEeBT4oZlZ7KulwB7A5sCJZvbJm6riRUd2P99D1vkws7skfYnwwoRNCReRR8T2DYrtf8DMzknU3Ts7FjO7TdLuse19CN+9PJ5wETQcuElSPbC3mdU38Ty4Fjry2sc26vHq6uqZUrNxj7mxdYUYoXCcla/V8VHDuvXK6tc2cO70Km5Z9FpbN69DSWMCvxW4ME6bDwWm8umL5S8A5prZifF1fouAB4ELgeFmdjqEKXBgV8IbkvoCNZKmxPpOAPYkJMeFkuYTksxRwDBCny0hRwKPdZ8GnEVITvsnypfGUXYuN0laA3wBONPMGuIU8CXA7oQXCNwvaWyMaYNyM5sg6fTMMSQdAfzLzA6Oy1uY2XvxAmWUmb0dj90beMrMLozbPWNmF8XPNwKHAHdntjWzr8W3K00FBmcCMLOPJGX38++yz4ekB4FTgSvN7KaY+LsR3gY1OE8ffTM7Fkk9CBcHh5nZW5KOJLw+8URJpwNnxxc6ZJ+fUwgXN/Tv35+Kioo8p6RzWLlyZbvEWFe3ca+ZGhoaqKur26jH3Ni6QoxQOM7s5J0s7wr9k5S6BG5mVZIGEkbf92atPhAYk7g/3ZP13/ObNDu+aWmNpBWE0elI4E4zWwUgaQbh4mCTWL46lud7Ew9m9n/A/0k6BvglcHwsL2skrGPN7ElJ/YFHJf2DMB1dEV/xh6SbgH0Jo/pc5TOz6qwGLpN0CXCPmT2S59gNhHf9ZoySdC5hmn0rwmg+k8BvibE8rPCaxX5mVtdIXPnOx2PABQrvHZ5hZi9IOV8tnDcWSYMJFxAPxH27Ed6k1Sgzuw64DqC0tNQ6+z9Xaa9/krOxD9kV/ulRV4gRCsc5YvJcanNcIA7oV8Kc8/bPsUfHowmFtylG2u6BZ8wivFz+lqxyAUeYWVn82SHxTttsaxKfGwgXM41lkaZOh98KjG3KDjEpL+HTGYBcGs10ibqeJ4zSq4FJcXScy4dmlpmy70l4ZeO3zWwI4ZWLPZPVZh+mQDNyng8zu5nwmsh6YE6ctm9qLAKeTtQ9xMxyvqPYOdd5nDO6lJIe6z+uU9KjG+eMLm2nFrWftCbwqcBFZladVT4HOENxSCZpWCz/gDBVXsjDwFhJveJ913HAI7F8nKQSSX0J94g3IOkLicWDgReKDSju34swTf8SsBDYT9I28eGyo4H5jZQDrI1Ty8Qp+NVm9jfCxc5X4jaN9UUmWb8tqQ/w7az1R8a6RwLvmdl7Weuz6855PiR9Hng5Pncwi3DrIm+78sRSA/SXtHfcpke8t14oRudcio0dNoBJhw9hQL8SRBh5Tzp8CGOHDWjvpm10qZtCBzCzN4Arc6y6GLgCqIpJYznhHu48YEJ8QGpSjv0y9S6RNI1wnxnCQ2yVAJJuIzzE9SohqedyuqQDgLWE+9PHZ1YUcQ+8HtgMmGZmi+M+58e2C7jXzO5qrJwwNVwlaQnwV+BSSetie36U2OY+SW+a2ais+Osk/Zkw0l0OPJHVznclPUp8iC1HHNn9nO98HAl8T9Ja4N+Ei7F3JC2Q9BRwX/IhNsJDcOvFEu+5fxu4StIWhN/lKwhT/tOAa/whNuc6p7HDBnTJhJ1NZk2dGXZdUXwKPeeDYWlVWlpqNTU17d2MNuX3TTuPrhAjdI04JS02s+EtrSetU+jOOedcl5bKKXS38ZlZeXu3wTnn3Kd8BO6cc86lkCdw55xzLoU8gTvnnHMp5AncOeecSyFP4M4551wKeQJ3zjnnUsgTuHPOOZdCnsCdc865FPIE7pxzzqWQJ3DnnHMuhTyBO+eccynkCdw555xLIU/gzjnnXAp5AnfOOedSyBO4c845l0KewJ1zzrkU6pAJXJJJujGx3F3SW5LuKbBfmaSDEssTJZ3dgnbk3F/SWZKekVQl6SFJOxaoZ5akp4o43kBJ9ZKWSlom6VFJpc1s+3p90cjxjmlO/S0laaikxyQ9LalaUs9Y/o8Y+9OSrpHULZY3qc+dc25mZS0jJs9lpwmzGTF5LjMra9u7Sa2qQyZwYBUwWFJJXP4GUEzPlwGNJq1WUgkMN7OhwHTg9/k2lHQ4sLIJdb9kZmVm9mXgBuAXzWxjGYX7YiCQM4FL6t7M4xYU6/4bcKqZfQkoB9bG1d+NsQ8G+gPfieUF+1zSeEkT26rdzrn0mFlZy/kzqqmtq8eA2rp6zp9R3amSeJv9kW4F9wEHE/5YHw3cAuwDIKk3cDUwhBDDxLj9RUCJpJHApFjPbpIqgB2AK8zsqljHWcCJcZvrzeyKWH4BcBzwOvAWsDi7YWY2L7H4OPC9XAFI6gOcBZwC3N608AHYHHg31tUTmAIMBz4GzjKzebnKgQVs2Bf/Bq7MhADsC0wGvihpKeFi4V1Cn/cEeksaA9wFbAn0AH5pZndJGgj8A1gIDAOeB44zs9VFxnUgUGVmywDM7L+ZFWb2fvzYHdg0trXoPnfOtZ4jr31sox+zrq6eKTUtP27la3V81LBuvbL6tQ2cO72KWxa91uL6O4KOnMBvBS6M0+ZDganEBA5cAMw1sxMl9QMWAQ8CFxJGaadDmAIHdgVGAX2BGklTYn0nAHsCAhZKmk+YkTiKkJS6A0vIkcCznES4eMjlYuB/gfUSW0yMw83swhz7DIoJtS/QK7YR4DQAMxsiaVfgfkm75CoHdsnRF3cDp5nZgnhh8SEwATjbzA6J24wH9gaGmtk7caQ8zszel7QN8LikWbE9pcBJsb6pwI+Bywr0VcYugEmaQxhl32pmn4yoY/kehH6dnmP/xvq8UZJOIVxQ0b9/fyoqKppTTWqsXLmy08cIXSPO9oixrq5+ox4PoKGhgbq6uhbXk528k+WtUX9H0GETuJlVxZHe0cC9WasPBMYk7k/3JIywc5ltZmuANZJWANsCI4E7zWwVgKQZhIuDTWL56lg+K0+dxPXfI4x898uxrgzY2cx+FuNIxjYLyFf3S2ZWFus4ErgO+GZs89Vx/+ckvUpIhPnKsy0ALpd0EzDDzN6QlOv4D5jZO5kwgN9J2hdYBwwg9B/A62a2IH7+G/ATik/g3WO7v0q4uHlI0mIzeyjGMTrOLNwE7A88kNkxu88lbQ08FFdvBWwqaWxc/r6ZVScPbGbXEfqU0tJSKy8vL7LJ6VRRUUFnjxG6RpztEWN7dGlrxTli8lxqc1yADOhXwpzz9m9x/S2hCa1TT0e9B54xi5AUbskqF3BEvFdcZmY7mNmzeepYk/jcQEgeOTNXZMU0TNIBhJmAMfECIdvewO6SlgP/BHaJU/lNMYsw1Q3529xYLJ8ws8nAyUAJYSS9a55NVyU+H0sYIe8eLyr+Q7hYgg37qah+i94A5pvZ2/Fi6V7gK1nt/ZAQ/2GZslx9bmb/zfweEGYdrkn8XqyXvJ1zXcc5o0sp6dFtvbKSHt04Z3SzngvukDp6Ap8KXJTjD/Ec4AzFIaSkYbH8A8LUcyEPA2Ml9Yr308cBj8TycZJKJPUFDs21czzetYREsiLXNmY2xcw+a2YDCaPN582svIi2JY0EXkq0+dh4/F0IMw41jZSv1xeSBplZtZldAjxJuLVQqL+2AFaY2VpJo4Dkk987SNo7fj6acJGCpEmSxhWIaw4wNPZ/d8Jo+hlJfSRtF+vpTngI77m4XLDPnXMuY+ywAUw6fAgD+pUgwsh70uFDGDtsQHs3rdV02Cl0ADN7g08fvEq6GLgCqIpJfDlwCDAPmBDvIU/KsV+m3iWSphHunUN4iK0SQNJtwFLgVUJSz+VSoA9wR7yGeM3MxsT9l2amwPMp8h64gI8Io2aAPwHXSKomPKw23szWSMpXnt0XI2MSbgCeIdxDXgd8LGkZMI34wFzCTcDdkp6MffJcYt2zwPGSrgVeIDxIB+HBwg1uDyRjNrN3JV0OPEEYud9rZrMlbQvMkrQZ0A2YC1wTq8jb5845l8vYYQM6VcLOJrOmzHw6F/79OHCPmQ3OsW6OmY3e+K1qutLSUqupqWnvZrSprnBvGLpGnF0hRugaccZnfoa3tJ6OPoXuUiYtyds559KuQ0+hu47JzJYTvmjFOedcO/ERuHPOOZdCnsCdc865FPIE7pxzzqWQJ3DnnHMuhTyBO+eccynkCdw555xLIU/gzjnnXAp5AnfOOedSyBO4c845l0KewJ1zzrkU8gTunHPOpZAncOeccy6FPIE755xzKeQJ3DnnnEshT+DOOedcCqUygUsySTcmlrtLekvSPQX2K5N0UGJ5oqSzW9COnPtLOkvSM5KqJD0kacc8+/9D0jJJT0u6RlK3AscbKKle0tK436OSSpvZ9vX6opHjHdOc+ltK0lBJj8W+qZbUM5bn7LNi+9w55zqLVCZwYBUwWFJJXP4GUFvEfmVAo0mrlVQCw81sKDAd+H2e7b5rZl8GBgP9ge8UUfdLZlYW97sB+EUz21hG4b4YCORM4JK6N/O4BcW6/wacamZfAsqBtXF1vj4rts+dc13AzMpaRkyey04TZjNi8lxmVhaTItIlrQkc4D7g4Pj5aOCWzApJvSVNlfSEpEpJh0naFLgIODKOYI+Mm+8mqULSy5J+kqjjLElPxZ8zE+UXSKqR9CCQc/RrZvPMbHVcfBzYPs9278eP3YFNAWtiH2wOvBvb1VPSX+JotVLSqHzlufpC0n7x89K4XV9gMrBPLPuZpPGS7pB0N3C/pD5xtLsk1n9YPOZASc9JuiGOiKdL6tWEuA4EqsxsWeyn/5pZQ2N9VmyfO+c6v5mVtZw/o5raunoMqK2r5/wZ1Z0uibfZKGojuBW4ME6bDwWmAvvEdRcAc83sREn9gEXAg8CFhFHa6RCmwIFdgVFAX6BG0pRY3wnAnoCAhZLmEy54jgKGEfpuCbC4QDtPIlxs5CRpDrBH3GZ6LBsT23lhjl0GSVoa29srthHgNAAzGyJpV0KC3SVXObBLjr64GzjNzBZI6gN8CEwAzjazQ+I244G9gaFm9k4cKY8zs/clbQM8LmlWbE8pcFKsbyrwY+CyAn2VsQtgsW/6A7ea2Scj6lx9lqXRPnfO5XfktY+16/Hr6uqZUtOyNlS+VsdHDevWK6tf28C506u4ZdFrLaq7I0ltAjezKkkDCaPve7NWHwiMSdyf7gnskKeq2Wa2BlgjaQWwLTASuNPMVgFImkG4ONgklq+O5bPy1Elc/z1gOLBfI3GMjvd3bwL2Bx4ws1lAvrpfMrOyWP+RwHXAN2Obr451PifpVUIizFeebQFwuaSbgBlm9oakXMd/wMzeyYQI/E7SvsA6YACh/wBeN7MF8fPfgJ9QfALvHtv9VWA18JCkxWb2UIxjgz7L7FiozyWdApwC0L9/fyoqKopsUjqtXLmy08cIXSPOjRVjXV19mx+jMQ0NDdTV1bWojuzknSxvad0dSWoTeDSLkBTKga0T5QKOMLOa5MaS9mRDaxKfGwh9kjNzRUVNc0s6gDATsF+8QMhfodmH8WLgMBLJqAizgL9kDpmvKcVUZGaTJc0m3Bd/PLY/l1WJz8cSRsi7m9laScsJF0uwYT815fbAG8B8M3sbQNK9wFeAhxLt3aDPiulzM7uOcNFDaWmplZeXN6FZ6VNRUUFnjxG6RpwbK8b27sbWiHPE5LnU5rgQGdCvhDnn7d+iuluDJrROPWm+Bw5h2vwiM6vOKp8DnKE4hJQ0LJZ/QJh6LuRhYKykXpJ6A+OAR2L5OEkl8R7xobl2jse7FhhjZivybNNH0nbxc3dC4nyuiLYljQReSrT52FjfLoQZh5pGytfrC0mDzKzazC4BniTcWijUX1sAK2LyHgUkn/zeQdLe8fPRwD/jcSZJGlcgrjnA0Nj/3Qmj6Wca67Ni+tw51zWcM7qUkh7r/6Oekh7dOGd0s/7RToeV6hG4mb0BXJlj1cXAFUBVTOLLgUOAecCEeA95UiP1LpE0jXDvHOB6M6sEkHQbsBR4lZDUc7kU6APcEa8hXjOzMXH/pXEKvDcwS9JmQDdgLnBN3KaYe+ACPgJOjuV/Aq6RVA18DIw3szWS8pVn98XImIQbgGcI95DXAR9LWgZMIz4wl3ATcLekJ2OfJC9AngWOl3Qt8AIwJZYPIcftgWTMZvaupMuBJwgj93vNbLakbfP1WWN97pzrWsYOGwDApXNq+FddPZ/tV8I5o0s/Ke8sZNbUB5+da1x8NuEeMxucY90cMxu98Vu1odLSUqupqSm8YYp1hall6BpxdoUYoWvEGZ/pGd7SetI+he5SpqMkb+ecS7tUT6G7jsnMlhO+aMU551wb8RG4c845l0KewJ1zzrkU8gTunHPOpZAncOeccy6FPIE755xzKeQJ3DnnnEshT+DOOedcCnkCd84551LIE7hzzjmXQp7AnXPOuRTyBO6cc86lkCdw55xzLoU8gTvnnHMp5AncOeecSyFP4M4551wKeQJ3zjnnUqhVErgkk3RjYrm7pLck3VNgvzJJByWWJ0o6uwXtyLm/pLMkPSOpStJDknbMsU0vSbMlPSfpaUmTizjeQEn1kpZKWibpUUmlzWz7en3RyPGOaU79aZXjd2SMpAnxc4t+X5xzLs1aawS+ChgsqSQufwOoLWK/MqDRpNVKKoHhZjYUmA78Ps92l5nZrsAwYISkbxVR90tmVmZmXwZuAH7RzDaWUbgvBgI5E7ik7s08bquT1K0Vqysj0S9mNsvMCl5cOec6tpmVtYyYPJedJsxmxOS5zKwsJmW4pNacQr8PODh+Phq4JbNCUm9JUyU9IalS0mGSNgUuAo6MI9gj4+a7SaqQ9LKknyTqOEvSU/HnzET5BZJqJD0I5Bz9mtk8M1sdFx8Hts+xzWozmxc/fwQsybVdAZsD78Z29ZT0F0nVMeZR+cpz9YWk/eLnpXG7vsBkYJ9Y9jNJ4yXdIelu4H5JfeIMw5JY/2HxmAPjzMINcRZiuqRejQUiaZqkayQ9Iul5SYfE8m6SLo3nskrSD2N5uaR5km4GquN2l8V2VEk6I263u6T5khZLmiNpu1heIekSSYvi8fbJ0y/jJf0xR3sHSfpHrPcRSbs28dw55zaSmZW1nD+jmtq6egyoravn/BnVnsSbqDVHbbcCFypMmw8FpgL7xHUXAHPN7ERJ/YBFwIPAhYSR8ekQpkSBXYFRQF+gRtKUWN8JwJ6AgIWS5hMuQI4ijJi7E5Lu4gLtPIlwsZFXbOOhwJVxeUxs54U5Nh8kaWlsb6/YRoDTAMxsSEwm90vaJVc5sEuOvrgbOM3MFkjqA3wITADONrNMMh0P7A0MNbN34ih8nJm9L2kb4HFJs2J7SoGTYn1TgR8DlxXoq4HAfsAgYJ6knYHjgPfM7KuSNgMWSLo/br8HMNjMXpH0I2AnYJiZfSxpK0k9gKuBw8zsrXjR9lvgxLh/dzPbI06Z/8rMDpCU3S/j87T1OuBUM3tB0p7An4D9C8TnXKOOvPax9m7CJ+rq6plS03Ha0xKVr9XxUcO69crq1zZw7vQqdtqcThNnW2u1BG5mVZIGEkbf92atPhAYo0/vV/YEdshT1WwzWwOskbQC2BYYCdxpZqsAJM0gXBxsEstXx/JZeeokrv8eMJyQlPJt050we3CVmb0cY5sF5Kv7JTMri/seSUgk34xtvjru/5ykVwmJOl95tgXA5ZJuAmaY2RuSch3/ATN7J9N84HeS9gXWAQMI/QfwupktiJ//BvyEwgn8djNbB7wg6WXCxdWBwFBJ347bbAF8AfgIWGRmr8TyA4BrzOzjGOs7kgYDg4EHYizdgDcTx5sR/7uYcPFQlHiB8zXgjkQfbZZn21OAUwD69+9PRUVFsYdJpZUrV3b6GKHt4qyrq2/1OpuroaGBurq69m5Gq8hO3snyhgbrNHG2tda+bzqLkBTKga0T5QKOMLOa5MZxpJRtTeJzQ2xjzswVWTENk3QAYSZgv3iBkM91wAtmdkUx9WaZBfwlc8h8TSmmIjObLGk24f7v47H9uaxKfD4W6A/sbmZrJS0nXCzBhv1UTL/l2kfAGWY2J7lCUnlWW5RjfwFPm9neeY6XOS+Z816sTYC6zIVUY8zsOsI5prS01MrLy5twmPSpqKigs8cIbRdnR+q6znQuR0yeS22Oi6MB/Ur45V6bdJo48wmP4bZca/8zsqnARWZWnVU+BzhDcXgkaVgs/4Aw9VzIw8BYhSfFewPjgEdi+ThJJfEe8aG5do7HuxYYY2Yr8h1E0m8II8ozi2hTLiOBlxJtPjbWuwthxqGmkfL1+kLSIDOrNrNLgCcJo99C/bUFsCIm71FA8mn7HSRlEufRwD/jcSZJGpenvu9I2kTSIODzsZ1zgB/F6XAk7RLPSbb7gVPjjAaStor798+0Q1IPSV9qJB6KiBkzex94RdJ3Yr2S9OUC9Trn2sk5o0sp6bH+s64lPbpxzuhm/SOeLqtVE7iZvWFmV+ZYdTHQA6iS9FRcBphHeGgt+RBbrnqXANMI984XAtebWWUsvw1YCvydkNRzuRToQ5hiXZqcao/3r5G0PWGEvhuwJG53clw3RtJFeeoeFLddBvwOODmW/wnoJqk6tnF8HPnnK8/uizMVHthbBtQT7ttXAR8r/JO1n+Voy03AcElPEi4SnkusexY4XlIVsBUwJZYPAf6dJ7YaYH489qlm9iFwPfBM7KOnCBdGuUbL1wOvEc75MuCY+HDgt4FLYtlSwtR3Y4r6HYnxnhTrfRo4rEC9zrl2MnbYACYdPoQB/UoQYeQ96fAhjB02oL2blioyK2oG2qVYfDbhHjMbnGPdHDMbnaN8Wtxnetu3sH2UlpZaTU1N4Q1TrDNNuzamK8TZFWKErhGnpMVmNryl9fg3sXVxuZK3c865jq/DfPmHaztmtpzw9HdT9hnfJo1xzjnXKnwE7pxzzqWQJ3DnnHMuhTyBO+eccynkCdw555xLIU/gzjnnXAp5AnfOOedSyBO4c845l0KewJ1zzrkU8gTunHPOpZAncOeccy6FPIE755xzKeQJ3DnnnEshT+DOOedcCnkCd84551LIE7hzzjmXQp7AnXPOuRQqmMAlmaQbE8vdJb0l6Z4C+5VJOiixPFHS2c1taL79JZ0l6RlJVZIekrRjnv1/K+l1SSuLPN5ASfWSlkpaJulRSaXNbPt6fdHI8Y5pTv3NaM++kpZI+ljStzfGMQuRVF7odyrHPrtIulfSi5KelXS7pG2bU5dzbuOaWVnLiMlz2WnCbEZMnsvMytr2blLqFDMCXwUMllQSl78BFNPTZUCjSauVVALDzWwoMB34fZ7t7gb2aGLdL5lZmZl9GbgB+EUz21hG4b4YCORM4JK6N/O4+bwGjAdubuV6NxpJPYHZwBQz29nMvghMAfq3b8ucc4XMrKzl/BnV1NbVY0BtXT3nz6j2JN5ExSaG+4CDCQnyaOAWYB8ASb2Bq4Ehsb6JcfuLgBJJI4FJsZ7dJFUAOwBXmNlVsY6zgBPjNteb2RWx/ALgOOB14C1gcXbDzGxeYvFx4Hu5AjCzx2OdRYa8gc2Bd2MdPQnJYjjwMXCWmc3LVQ4sYMO++DdwZaZpwL7AZOCLkpYSLhbeJfR5T6C3pDHAXcCWQA/gl2Z2l6SBwD+AhcAw4HngODNbnS8QM1se41jXWMCSvgP8CmgA3jOzfePxbgR6x81ON7NHJZUDvwb+Q7hgmQFUAz8FSoCxZvaSpGnAh8CXgG1j3603Ws71O2Vmd2U17xjgMTO7OxHXvLh/eWNxOddWjrz2sTaru66unik1bVf/xlT5Wh0fNaz/56d+bQPnTq9ip83pNHG2tWIT+K3AhXFacigwlZjAgQuAuWZ2oqR+wCLgQeBCwsj4dAhT4MCuwCigL1AjaUqs7wRgT0DAQknzCbMDRxGSUndgCTkSeJaTCBcPRYuJcbiZXZhj9aCYUPsCvWIbAU4DMLMhknYF7pe0S65yYJccfXE3cJqZLZDUh5DQJgBnm9khcZvxwN7AUDN7J47Cx5nZ+5K2AR6XNCu2pxQ4KdY3FfgxcFlT+iGPC4HRZlYbzy3ACuAbZvahpC8QLuaGx3VfBr4IvAO8TLgY20PST4EzgDPjdgOB/YBBwDxJO2cdd4PfKUkPmtmqxDaDKfz7sAFJpwCnAPTv35+KioqmVpEqK1eu7PQxQseJs66uvs3qbmhooK6urs3q35iyk3eyvKHBOk2cba2oBG5mVXHkdTRwb9bqA4ExifvTPQkj7Fxmm9kaYI2kFYQR2EjgzswfZ0kzCBcHm8Ty1bF8Vp46ieu/R0gk+xUTUyK2WUC+ul8ys7JY/5HAdcA3Y5uvjvs/J+lVQqLOV55tAXC5pJuAGWb2Rp6ZgQfM7J1MiMDvJO0LrAMGEPoP4HUzWxA//w34Ca2TwBcA0yTdThhRQxj9/1FSGWFknozvCTN7E0DSS4QLGAgj8VGJ7W43s3XAC5JeJlzYJeX7nXq2pQGZ2XWE80hpaamVl5e3tMoOraKigs4eI3ScONuyCR0lxtYwYvJcanNc7AzoV8Iv99qk08SZjya0Tj1NeQp9FiEp3JLdFuCIeK+4zMx2MLN8f2jXJD43EC4gGpvTtmIaJukAwqhtTLxAaAuzCFPdkL/NRc3Pm9lk4GTC1PLjcbSeS3LEeSzh/u7u8aLiP4TEBhv2U1H9li0+6Lc0zjpgZqcCvwQ+ByyVtDXws3jsLxMumDZNVJHs+3WJ5XWsf7FYqL3F/E49DezelPiccx3DOaNLKenRbb2ykh7dOGd0s54T7rKaksCnAheZWXVW+RzgDMUhpKRhsfwDwtRzIQ8DYyX1ivc+xwGPxPJxkkok9QUOzbVzPN61hOS9ognxNNVI4KVEm4+Nx9+FMDqsaaR8vb6QNMjMqs3sEuBJwgi0UH9tAawws7WSRgHJp+13kLR3/Hw08M94nEmSxhUboJldkEmaiXYujLcX3iYk8i2AN+MI+vtAt7wV5vcdSZtIGgR8ntBHSfl+p5JuBr4m6eBMgaRvShrSjPY45zaiscMGMOnwIQzoV4III+9Jhw9h7LAB7d20VCn66WYze4NPH7xKuhi4AqiKf3CXA4cA84AJcTQ3Kcd+mXqXxAebFsWi682sEkDSbcBS4FVCUs/lUqAPcEf8e/+amY2J+y9NJKPfEx586iXpjXiciUXeAxfwEWHUDPAn4BpJ1YSH1cab2RpJ+cqz+2JkTMINwDOE+/brgI8lLQOmER+YS7gJuFvSk7FPnkusexY4XtK1wAuEB+kgPAS2we0BSV8F7iQ8EHeopF+b2Zdy9W28zy3gIWBZjP3v8QG3eaw/S1CsGmA+4RbAqfF+enJ9vt+pT5hZvaRDgCskXQGsBaoID81t3Yw2Oec2orHDBnjCbiGZNWu21XUQ8dmEe8xscI51c8xs9MZvVX7xYu0eM5ve3m0pLS21mprswX/n0pnumzamK8TZFWKErhGnpMVmNrzwlo3zb2LrxDpa8nbOOdd6WvsLQtxGFv9N9waj747KzMa3dxucc64z8BG4c845l0KewJ1zzrkU8gTunHPOpZAncOeccy6FPIE755xzKeQJ3DnnnEshT+DOOedcCnkCd84551LIE7hzzjmXQp7AnXPOuRTyBO6cc86lkCdw55xzLoU8gTvnnHMp5AncOeecSyFP4M4551wKeQJ3zjnnUqjVE7gkk3RjYrm7pLck3VNgvzJJByWWJ0o6uwXtyLm/pH0lLZH0saRvF1lXhaQnE8vDJVU0t23N0dL+aGuSxkia0Ep1nSmpV2L5Xkn9WqNu51z6zKysZcTkuew0YTYjJs9lZmVtezepQ2iLEfgqYLCkkrj8DaCY3i4DDiq0USt4DRgP3NzE/T4j6VvNOaCk7s3Zrz1J6taU7c1slplNbqXDnwl8ksDN7CAzq2ulup1zKTKzspbzZ1RTW1ePAbV19Zw/o9qTONBWieU+4GBgOnA0cAuwD4Ck3sDVwJB4/Ilx+4uAEkkjgUmxnt3iSHcH4AozuyrWcRZwYtzmejO7IpZfABwHvA68BSzObpiZLY/brmtiTJcCv4xt/YSknsAUYDjwMXCWmc2TND72QU+gt6S/AmOBbsBg4H+BTYHvA2uAg8zsHUk/AE6J614Evm9mq/M1StKhsV2bAv8FjjWz/0iaCAwCBgCfA35vZn+WVE7o6/8CpcDDwI/NbJ2klcDlwGjg55L2IKufY98PNrMTJQ0hnNs9gO8Cw83sdEnTgHpgV2BH4ATgeGBvYKGZjY9tnwJ8FSgBppvZryT9BPgsME/S22Y2StLyWPfbuc69pIHxvPwT+BrhgvEwM6vP12/ONdeR1z7WLsetq6tnSk37HHtjyo6z8rU6PmpY/891/doGzp1exS2LXtvYzetQ2iqB3wpcGKfNhwJTiQkcuACYGxNAP2AR8CBwITEBQJgyJiSAUUBfoCb+wR9KSAh7AgIWSppPmE04ChgW41pCjgTeGElLzawsz+rHgHGSRgEfJMpPAzCzIZJ2Be6XtEtctzcwNCbm8YTEPYyQ1F8EzjOzYZL+QLjwuAKYYWZ/ju35DXAS4YInn38Ce5mZSToZOBf4eVw3FNgL6A1USpody/cAdgNeBf4BHE642OoNPGVmF0randz9fAVQIWkc4Vz+0MxWS8pu15bA/sAY4G5gBHAy8ISkMjNbClwQ+6Yb8JCkoWZ2VUzSo8zs7WSFjbTpXeALwNFm9gNJtwNHAH/LbpSkUwgXSPTv35+KiopGujb9Vq5c2eljhI0bZ11d+1wXNjQ0UFdX1y7H3piy48xO3snyrtAfjWmTBG5mVXFUdDRwb9bqA4Exifu5PQkj7Fxmm9kaYI2kFcC2wEjgTjNbBSBpBuHiYJNYvjqWz2pGu8sKbPIbwmj3vETZSGKCNbPnJL0KZBL4A2b2TmLbeWb2AfCBpPcIiQ2gmpBsIdx++A3QD+gDzCnQpu2B2yRtRxiFv5JYd1cchdZLmkdI3HXAIjN7GUDSLTGG6UAD8PdEXBv0s5lVxouRKuBaM1uQp113x4uKauA/ZlYd63kaGAgsBb4bE2p3YDvCRUVVI7HmO/ezgFfiRQGEC7eBuSows+uA6wBKS0utvLy8kcOlX0VFBZ09Rti4cbZXd3bVczli8lxqc1w0DehXwpzz9t+ILWs9rfO0UNs+hT4LuIwwxZok4AgzK4s/O5jZs3nqWJP43ED4Q7/BUC/Bmt3aIpjZXMIFx16J4sbasyprORnPusTyOj69mJoGnG5mQ4Bfx+M15mrgj3H7H2Ztn90fVqD8QzNriJ8bi+sLwErCVHc+ydiy4+4uaSfgbODrZjYUmE3hWBtrU67fFedcyp0zupSSHus/klPSoxvnjC5tpxZ1HG2ZwKcCF2VGXglzgDMU51wlDYvlHxCmygt5GBgrqVe8nz4OeCSWj5NUIqkvcGhrBJHDbwnT1Mn2HAsQp853AGpaUH9f4E1JPTL1FrAFnz4keHzWusMk9ZS0NVAOPBHL95C0k6RNgCMJ0/DZcvazpC2AK4F9ga2LfZI/h80JFzjvSdoWSD4gmO93Id+5d851UmOHDWDS4UMY0K8EEUbekw4fwthhA9q7ae2uzUYpZvYG4Q99tosJ91GrYhJfDhwCzAMmSFrKpw+x5ap3SXxIalEsut7MKgEk3UaYmn2VPH/YJX0VuJNwj/ZQSb82sy/FdY3dA88c/15JbyWK/gRcE6eKPwbGm9maHPeEi/U/wMIYQzWFL2omAndIqgUeB3ZKrFtEGNnuAFxsZv+KFxmPAZMJDxI+TOiP9eTrZ0lTgT+Z2fOSTiI8bPZwU4M0s2WSKoGngZeB5FT8dcB9kt40s1FFtGlgU4/vnEuPscMGeMLOQWZtOuvs2kl8CHClmV2WVV4OnG1mh7RDszqU0tJSq6lpyWRJx9dV75t2Rl0hRugacUpabGbDW1qPfxObc845l0I+AnddlqQPaNnzCmmwDfB2wa3SryvE2RVihK4RZ6mZFfPMV6P8SV3XldW0xjRWRybpyc4eI3SNOLtCjNA14lTiq7lbwqfQnXPOuRTyBO6cc86lkCdw15Vd194N2Ai6QozQNeLsCjFC14izVWL0h9icc865FPIRuHPOOZdCnsCdc865FPIE7jodSd+UVCPpRWnD9/5IOkfS0vjzlKQGSVvFdcslVcd1rfJPPdpKEXFuIeluScskPS3phGL37ShaGGNnOpdbSrpTUpWkRZIGF7tvR9HCGFNxLiVNlbRC0lN51kvSVbEPqiR9JbGu6efRzPzHfzrND9ANeAn4POH1qsuA3RrZ/lDC++kzy8uBbdo7jtaIE/gFcEn83B94J27bpD5KY4yd8FxeCvwqft4VeKjYfTvCT0tiTNm53Bf4CvBUnvUHAfcR3qy4F7CwJefRR+Cus9kDeNHMXjazj4BbgcMa2f5oNnzlbRoUE6cBfeNLg/oQktvHRe7bEbQkxjQpJs7dgIcAzOw5YGB8i19nOpf5YkwNM3uY8DuYz2HAXy14HOgnaTuaeR49gbvOZgDwemL5jVi2AUm9gG8Cf08UG3C/pMWSTmmzVrZcMXH+Efgi8C/Cm+1+ambrity3I2hJjNC5zuUy4HAASXsAOwLbF7lvR9CSGCE957KQfP3QrPPoX6XqOptc73HN928lDwUWmFnyinmEhdeufgZ4QNJz8aq6oykmztGE1+vuDwwixPNIkft2BM2O0czep3Ody8nAlQqvW64GKgkzDZ3pXOaLEdJzLgvJ1w/NOo8+AnedzRvA5xLL2xNGZ7kcRdb0uZn9K/53BeE96Xu0QRtbQzFxngDMiNN1LwKvEO4tNqWP2lNLYuxU59LM3jezE8ysDDiOcL//lWL27SBaEmOazmUh+fqhWefRE7jrbJ4AviBpJ0mbEpL0rOyNJG0B7AfclSjrLalv5jNwIJDzadIOoJg4XwO+DhDvJZYCLxe5b0fQ7Bg727mU1C+uAzgZeDjOMnSac5kvxpSdy0JmAcfFp9H3At4zszdp5nn0KXTXqZjZx5JOB+YQnuycamZPSzo1rr8mbjoOuN/MViV23xa4MzwPRXfgZjP7x8ZrffGKjPNiYJqkasIU3Xlm9jZArn3bI47GtCRGSZ+nc53LLwJ/ldQAPAOc1Ni+7RFHY1oSIyn6/1LSLUA5sI2kN4BfAT3gkxjvJTyJ/iKwmjCD1Ozz6F+l6pxzzqWQT6E755xzKeQJ3DnnnEshT+DOOedcCnkCd84551LIE7hzzjmXQp7AnXPtQuEtcJk3wt0Rv9q2uXVNk/Tt+Pl6Sbs1sm25pK8llk+VdFxzj+1ce/EE7pxrL/VmVmZmg4GPgFOTKyV1a06lZnaymT3TyCblwCcJ3MyuMbO/NudY7UWSf4eH8wTunOsQHgF2jqPjeZJuBqoldZN0qaQn4vuTfwifvFf5j5KekTQb+EymIkkVkobHz9+UtEThfeEPSRpIuFD4WRz97yNpoqSz4/Zlkh6Px7pT0paJOi9ReE/185L2yQ5AUp94jCUK764+LLHuuFjnMkk3xrJt4zGWxZ+vSRqoxLukJZ0taWKiDb+TNB/4qaRDJS2UVCnpQcU3d8V2/CW2oUrSEZJOkvSHRL0/kHR565w61178Ks45167iaPJbQObbtfYABpvZKwpvnnrPzL4qaTNggaT7gWGEr00dQvimrmeAqVn19gf+DOwb69rKzN6RdA2w0swui9t9PbHbX4EzzGy+pIsI36R1ZlzX3cz2kHRQLD8gK5QPgXHx6z+3AR6XNIvwmswLCC/keFvSVnH7q4D5ZjYuzjb0AbYs0F39zGy/2O4tgb3MzCSdDJwL/Bz4n9hnQxLbfQRUSTrXzNYSvgHshwWO5To4T+DOufZSovDmKQgj8P9HmNpeZGavxPIDgaGZ+9vAFsAXgH2BW8ysAfiXpLk56t+L8H3amRdiNPae5sz34/czs/mx6AbgjsQmM+J/FwMDc1UB/E7SvkDmta3bEt6UNj3zNbaJduxPeGkHMY73MiP+RtyW+Lw9cJvC+6Q3Jb74g3BhcVRmIzN7N8Y3FzhE0rNADzOrLnAs18F5AnfOtZf6+OapT8Tvu05+P70II+I5WdsdROHXLaqIbZpiTfxvA7n/dh5LeIPW7ma2VtJyoGcT2/Ex69/a7Jm1Ptk3VwOXm9ksSeXAxFie73jXA78AngP+UmR7XAfm98Cdcx3ZHOBHknoASNpF4Y1UDwNHxXvk2wGjcuz7GLCfpJ3ivpmp6w+Avtkbm9l7wLuJ+9vfB+Znb9eILYAVMXmPAnaM5Q8B35W0dVY7HgJ+FMu6Sdoc+A/wGUlbx1sGhxQ4Xm38fHyi/H7g9MxCZlRvZgsJr6w8hqzX6Lp08gTunOvIrifc314SH+66ljD6vRN4AagGppAj0ZrZW8ApwAxJy/h0+vluYFzmIbas3Y4HLpVUBZQBFzWhrTcBwyU9SRiNPxfb8TTwW2B+bEfm4bGfAqMU3qS2GPhSvD99EbAQuCdTRx4TgTskPQK8nSj/DbClwj/PW8b6Fze3Awsy0+ou3fxtZM4510VIugf4g5k91N5tcS3nI3DnnOvkJPWT9DzhuQNP3p2Ej8Cdc865FPIRuHPOOZdCnsCdc865FPIE7pxzzqWQJ3DnnHMuhTyBO+eccyn0/wGlgeUtk1AfqQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "labels = list(results.keys())\n",
    "\n",
    "means = np.array([results[k][\"Test accuracy\"] for k in labels])\n",
    "lower_error = np.array([results[k][\"Lower 95% CI\"] for k in labels])\n",
    "upper_error = np.array([results[k][\"Upper 95% CI\"] for k in labels])\n",
    "\n",
    "asymmetric_error = [means - lower_error, upper_error - means]\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(7, 3))\n",
    "ax.errorbar(means, np.arange(len(means)), xerr=asymmetric_error, fmt=\"o\")\n",
    "ax.set_xlim([0.75, 1.0])\n",
    "ax.set_yticks(np.arange(len(means)))\n",
    "ax.set_yticklabels(labels)\n",
    "ax.set_xlabel(\"Prediction accuracy\")\n",
    "ax.set_title(\"95% confidence intervals\")\n",
    "\n",
    "plt.grid()\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"matplotlib-figures/comparison.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- We can see that bootstrapping the test set results in the widest confidence intervals, and the .632 estimates results in the smallest confidence intervals. If the .632 confidence intervals are correct (contain the true parameter 95% of the time), then these would be most desirable. \n",
    "- However, which confidence interval method is most indeed correct is tricky to answer. There is a quick attempt in the notebook linked in the next section.\n",
    "- Note that the normal approximation and the bootstrap test set methods have slightly different test accuracy means (dots) because the bootstrap test set one is from averaging the bootstrap test set accuracies; however, we could use the original test set accuracy as the center."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence Intervals and the True Model Performance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Have a look at the supplementary [./ci-simulation.ipynb](./ci-simulation.ipynb) notebook for a quick simulation study."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus: Creating Confidence Intervals with TorchMetrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- In the context of deep learning and PyTorch, I recently wrote about [TorchMetrics](https://sebastianraschka.com/blog/2022/torchmetrics.html), which is an awesome tool for evaluating models in cases where the dataset is too big to fit into memory.\n",
    "- We can also use TorchMetrics for bootstrapping. In this case, bootstrapping the test set.\n",
    "- For comparison, the code for bootstrapping the test set we used earlier was as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.956304347826087"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.fit(X_train, y_train)\n",
    "\n",
    "predictions_test = clf.predict(X_test)\n",
    "acc_test = np.mean(predictions_test == y_test)\n",
    "\n",
    "rng = np.random.RandomState(seed=123)\n",
    "idx = np.arange(y_test.shape[0])\n",
    "\n",
    "test_accuracies = []\n",
    "\n",
    "for i in range(200):\n",
    "\n",
    "    pred_idx = rng.choice(idx, size=idx.shape[0], replace=True)\n",
    "    acc_test_boot = np.mean(predictions_test[pred_idx] == y_test[pred_idx])\n",
    "    test_accuracies.append(acc_test_boot)\n",
    "\n",
    "bootstrap_train_mean = np.mean(test_accuracies)\n",
    "bootstrap_train_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8695652173913043 1.0\n"
     ]
    }
   ],
   "source": [
    "ci_lower = np.percentile(test_accuracies, 2.5)\n",
    "ci_upper = np.percentile(test_accuracies, 97.5)\n",
    "\n",
    "print(ci_lower, ci_upper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Using the [`Bootstrapper`](https://torchmetrics.readthedocs.io/en/stable/references/wrappers.html#bootstrapper from TorchMetrics, we can we can replicate the results from above as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'mean': tensor(0.9602), 'std': tensor(0.0408), 'quantile': tensor([0.8696, 1.0000])}\n"
     ]
    }
   ],
   "source": [
    "import torch\n",
    "from torchmetrics import Accuracy, BootStrapper\n",
    "\n",
    "torch.manual_seed(123)\n",
    "\n",
    "quantiles = torch.tensor([0.05, 0.95])\n",
    "base_metric = Accuracy()\n",
    "bootstrap = BootStrapper(\n",
    "    base_metric, num_bootstraps=200, sampling_strategy=\"multinomial\", quantile=quantiles\n",
    ")\n",
    "\n",
    "bootstrap.update(torch.from_numpy(predictions_test), torch.from_numpy(y_test))\n",
    "output = bootstrap.compute()\n",
    "\n",
    "print(output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Note that in practice this can come really handy if we want to compute things incrementally (e.g., if the test set is too big for memory).\n",
    "- E.g., let's assume the predictions come in multiple chunks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = np.arange(predictions_test.shape[0])\n",
    "groups = np.array_split(idx, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'mean': tensor(0.9550), 'std': tensor(0.0421), 'quantile': tensor([0.8696, 1.0000])}\n"
     ]
    }
   ],
   "source": [
    "torch.manual_seed(123)\n",
    "\n",
    "quantiles = torch.tensor([0.05, 0.95])\n",
    "base_metric = Accuracy()\n",
    "bootstrap = BootStrapper(\n",
    "    base_metric, num_bootstraps=200, sampling_strategy=\"multinomial\", quantile=quantiles\n",
    ")\n",
    "\n",
    "for group in groups:\n",
    "\n",
    "    pred_chunk = torch.from_numpy(predictions_test[group])\n",
    "    label_chunk = torch.from_numpy(y_test[group])\n",
    "\n",
    "    bootstrap.update(pred_chunk, label_chunk)\n",
    "\n",
    "output = bootstrap.compute()\n",
    "print(output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
